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A General Method for the Computation of Cartesian Coordinates 


And Partial Derivatives of the Two-Body Problem 


A general solution for cartesian coordinates and partial derivatives 
of the two-body problem has been programmed as a double-precision 
Fortran 4 program for the IBM 7094. The compact subroutine provides 
an accurate and efficient computation of coordinates and partial derivative 
for all cases of two-body motion. A derivation of all equations used by 
the subroutine is given for those interested in the formulation. A de¬ 
scription of the subroutine is also given for those interested in the 
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INTRODUCTION 


The classical two-body problem may be reduced to the solution 
of the differential equations 


dr 

dt 


= r 


dr 

dT = 


where is a known constant and 


r= \J~fT* 

in the magnitude^ of r. A solution of the two-body problem gives the 
position vector r^ and velocity vector i^at any time t in terms of the 
position vector r^ and velocity vector r^ at a reference time t . Many 
practical applications also require the partial derivatives of the com- 

z of r and x. v. z of r with respect to the components 

r_. Classical formulas for the coordi- 


ponents x, y 

V V z o of ~o. ln ? 

nates x, y, z, x, y 


r and x, 

i • 

, y 


§’ 


qI z^ of 
and partial derivatives 


9 (x, y, z, x, y, z) _ 

8 <WVVV i o' 

have the disadvantage that several different formulations must be avail¬ 
able and one must be selected depending on the particular values of x 


O’ 


V V V z o ana " • °’ 

General formulas are available for computing the coordinates x, y,z, 


x, y,z and partial derivatives 


d (x, y, z, x, y, z) 

9 (x . y z„, x . 

0 0 0 0 


as well as 


o' z o> 


(x, y, z, X, y, z) r „ ... 

Tm - f ° r a11 P° ssible values of x Ql y , z Q , x , y , z t t and |i . 

The formulas have been programmed as a double precision Fortran 4 sub¬ 
routine for the IBM 7094. The subroutine is advantageous in that it offers 
an efficient and compact program for the accurate computation of coordinates 
and partial derivatives for all possible cases of two-body motion. The 
purpose of this report is to give a derivation of all equations used by the 
subroutine and a detailed description of the computations. For those who 
are not interested in these details but only in using the subroutine, the 
inputs and outputs are described in Paragraph 4. 1 below and a Fortran 4 
listing is given in Paragraph 4. 8. 


The solution for the coordinates is derived in Section 1 by defining a 
new independent variable ip by the differential equation 


dip 

dt 
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which Sundman applied to a theoretical investigation of the three-body 
problem. The solution is a modification of that given originally by 
Stumpff and differs primarily in that it is valid for all values of the 
constant fi . The relation of the general solution for coordinates to the 
classical formulas for elliptic two-body motion is also considered. It 
is shown that ip is a generalization of the eccentric anomaly and that 
the relation between t and <p is a generalization of Kepler's equation. 

The partial derivatives are obtained in Section Z by differentiating 
the solution for coordinates and manipulating the results algebraically. 
The approach is an extension of Sconzo's derivatives which are based 
on those of KUhnert. Although the derivation is tedious and involved, 
the final derivatives are compact and general. Some formulas used in 
the derivation of the coordinates and partial derivatives are discussed in 
Section 3 with formulas that are used for computation. A detailed 
description of the double-precision Fortran 4 subroutine for the IBM 7094 
is given in Section 4. The description should be sufficient for those 
wishing to make modifications or additions to the subroutine. 
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1. THE SOLUTION FOR THE COORDINATES 


1. 1 THE EQUATIONS OF MOTION OF THE TWO-BODY PROBLEM 

The equations of motion of the two-body problem are 

r. - r 


m r = -G m, rn 

11 12 


1 


12 


m r 
2 2 


- G n^m 


r 2 ' r i 


12 


where 


4 


12 


\fi*: i" 


V ' - ? 2 ) 


is the distance between the two bodies and G is the universal 
constant of gravitation. At the time t the body of mass m has the in- 

. i - 1 

ertia 1 cartesian position vector r^ and velocity vector ? , and 

the body of mass m has the inertial cartesian position vector r 

^ 2 

and velocity vector 7 Each body is attracted toward the other 

2 

with a force of magnitude G m^m^/r^ where the acceleration of 

_> s i 

m toward m is r and the acceleration of m toward m is 7 . 

1 £ 1 2 12 

Initial conditions are specified by the position 7 and velocity 

• 1,0 

1,0 of m and the position r and velocity 7 of m at a 

1 2,0 y 2,0 2 

given reference time t^. 

A closed-form solution of the differential equations above gives 
r^, and r r^ at any time t in terms of 7^ ^, 7^ ^ and 

r 2 O' r 2 0 t Q as we ^ as the constants m^, m^ and G. Such 

a solution is extremely important, for many practical applications 
in astronomy and astronautics. However, it is convenient to trans¬ 
form the differential equations by defining the position R of the 

center of mass and the position r of m relative to m by 

2 1 


1 



R = 


m i r i + m 2 r 2 
m i + m 2 


r = r — r 
2 1 


so that 


R = 


m l r i tm 2 r 2 
m l + m 2 


r = 


r — r 
2 1 


give the velocity R of the center of mass and the velocity r* of 

m relative to m,. Thus the values of the four quantities above 
2 1 

at the reference time t^ are 


R o = 


m, r, „ + m_ r 

1 1,0 _ 2 2,0 

m, + rru 

1 2 


V 


m l r i,0 + m 2 r 2,0 

m i + m 2 


and 


r = r — r 

0 2,0 1,0 


r = r — r 

0 2,0 1 , 0 . 


Time differentiation of R and r above and substitution from the 

differential equations show that 

• ♦ 

R = 0 


where 


"r = - G (m j + m^) r/ r 


r = VTT 


An integration of these differential equations gives R, R and r, 
at any time t. Simultaneous solution of the defining equations fo 
R and r shows that 


m 


? = R-- 

1 m. 4 m, 


m 


1 


1 


r = R 4 — , 

2 m i 4* m 


2 



and time differentiation gives 


m 


r i = R 


m i + m 2 


» 

r 


m 


i 


r~ = R + — 

2 m, + m. 

1 i 


r 


so that r , r and r , r may always be determined from R, R 

XX L* L* 

and r, r . 

Thus the integration of the differential equations for ? and r^ 

has been transformed to an integration of differential equations for 

fl and lr. Integration of the differential equation 
»• 

R = 0 

*■ > 

for R shows that 

R = t 

0 


and 


R = 


R o + 


R o “-‘o' 


so that the center of mass moves along a straight line at constant 

—* -V —► -U 

velocity. However, the inertial coordinates R^, R^ and thus R, R 
at any time t are usually unknown or the center of mass is chosen 
as the origin of inertial coordinates. In any case, the two-body 
problem reduces to the solution of the differential equation 

— 3 

r = - ii r / r 

where 

1 1 

r = v r * r 

is the magnitude of r and 
M = G (m^ + m^) 

is a postive constant. This differential equation is actually ob¬ 
tained for all inverse-square Newtonian forces between two bodies, 


3 



but the constant \i may be different. For example, the same 
differential equation describes the electrostatic repulsion of 
two electrons but ^ is then negative and is a more complicated 
function of the electronic mass and charge. 
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1.2 DERIVATIVES WITH RESPECT TO A NEW VARIABLE 


The differential equation 

—. 3 

r = - M r / r 


where 


r = r • r 

is conveniently integrated by d 
4 by the differential equation 
4) - 1/ r 

where 4 1 is zero when t is t^. 
with respect to ip , 


efining a new independent variable 
of the Sundman transformation, 

Using a prime to denote differentiation 


t' = 1/4 = 1 / (1 / r) 


= r 

from the properties of derivatives. The chain rule for differentiation 
shows that 


= r r. 

The method of solution described 
of the successive derivatives of t and 
second derivative of t is 


below requires 
r with respect 


determination 
to ^ . The 


t 




/ 

r 


2 


( 7 * 7 ) 2 





r 


• r 


where the quantity 
a = "f • 

is defined so that 

t" = r / = a . 

The second derivative of r*is 

-+// -V ' . . -V , 

r = r r + r r 


5 



but 


"r^ = ? t / = (- r 

= - M r/ r 

so that 

r = — /x r / r + ^r. 

The third derivative of t is 

. /// // _ / -*7 -V . -► -^/ 

t = r = <r = r • r + r*r 

• r * r t t ? 

= rr • r + r • (- n r/r ) 

9 

= r (r • ~r) - I* . 

However, the quantity 

Of 3 r • r - 2 /i/r 
is defined to obtain 

t'" = r H = a /= r(T •’r) -2 n + fi = rfr • 7 -2 /x /r) 

= a r + /i . 

The third derivative of r is then 

-*/ . ^ 2 /-v -v / 

= - ji r /r + ptrr/r + cr r + ar 

-V -► 2 -V 2 

= - /i r + /n r <r/r + (a'r+n)r-cr,jxr/r 


+ M 


r 


= a r r. 

The quanity a has the important property that 

• * 9 

c/ = 2r • r* + 2 M r'/r 

-*► 2 2 
= 2r * ( - M r/r ) + 2 |i cr/r 

= -2 M a/r^ + 2/i cr/r 2 


= 0 


so that C* is constant for every value of $ as well as t. This 
property gives all the derivatives of t and r as simple functions 
of the first two derivatives t 1 , t H and r’, r". From the results 

above, 

t m = Q? t' + M 
= a ? 


6 



so that 


t ,m = a t" 


r HH — Q| J. 11 


JTtim 


t ,,mi 

r nun 
^ n inn 


: a t m = a ( ot t* + M ) 

2 

a t' + a n 

a r" 1 = a ( a r^) 

2 

Q! r 1 
2 

a t" 

2 

a r 11 

a ^ t m = a ^ ( a t* + n ) 

3 2 

= Of t 1 + Of /i 


2 

a r 1,1 

3 - 
a r ? 


o: (Of r 1 ) 


etc. 

The final result of the above derivation is that 
t' = r 
t" = cr 


111 - 

Of 

r 

+ M 

till - 

Of 

a 


urn - 

a 

2 

r 

+ a: 

niiiL- 

ot 

2 

c 

F 

««tt in 

5 

= a 

r + 

tninn 

= a 3 

(j 


etc. 


7 



and that 


r 1 

= 

♦ 

r r 

7 " 


• 

- n r / r + a 7 


- 

a ( r r) 

■fim 

= 

♦ 

O' (- (i "r / r + a "r ) 


= 

2 , ^ 

O' (r r) 

r 1 lift 1 

- 

Z < -* v 

-O' (- M r/ r + cr r) 

r mtm 

- 

O' (r r) 


- 

3 , * 

Of (- /i r/r + a r) 


etc. 


where 


a - r 


/ . 


0 


0 


-2 [it 


0 


since ot is zero and a is invariant for all values of 4 as well as t. 
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1.3 SERIES SOLUTION IN THE NEW VARIABLE 


The Taylor's series for t in terms of >P is 

t = t Q + t 0 ' + t Q " tp 2 /2 : + t 0 "' / 3 ! + t "" 4> 4 /4l . 

where t Q , t Q \ tg", tp"', t "". are the values of t, t\ t", 

t" 1 , t"", ... when i p is zero, i. e. when t is t . But 


t 1 
0 


= r 


0 


t ” = a 

0 0 


V 


= a r 0 + n 


f nit 

o 


f mti 

0 


t """ 

0 


= a (t 


0 

Z 


O' 


r 0 + a M 


z 

O' a 


0 


etc. 

from the results above where 

“ = W 2 " /r o. 

Thus substitution for t ' t ", t m , t ,m , ... in the Taylor’s 

0 0 0 0 

series gives 


t = t + r $ 

0 o q 


+ CT 0 tp / 2 ! 


+ ( a r Q + n ) ip /3\ + a a 


, 4 , 

4> /4! 


+ ( a r ) ip /5! + a a ip /6! 


0 


0 


• • 


9 



or, collecting coefficients of r , cr^, and p , 


t = t 


0 


3 2 

+ i*q{ ^+a^/3!+a 


^ / 5 ! 4 a ^ ^ / 7 ! 4 . . 


+ CT ( ^ 2 /2 ! + a ip 4/41 + a 2 </> /6! + a 3 ^ 8 /8! + 


+ fi ( 4> 3 /3 ! + a i/ 5 /5! + 


7 3 9 

^/7!+ a $ /9! 4 


gives the value of t corresponding to any particular value of 


Similarly, the Taylor's series for r in terms of ip is 


where 


"r = T q + + 7" ip / 2! + r^'" 


^ / 3 ! + r "" 


ip /4! 4 . . 


• • • • 


0 

V 


= r r 
0 0 


= - (* V r 0 * *0 *0 


III 


0 


= a r. r 
0 0 


llll 


0 


“ (- M + "(jV 


0 0 


nm — 


0 


Of r 


0 0 


mm _ 


0 


“ "VV Vo» 


etc. 


from the results above. Thus substitution for "r \ "r ", r 
in the Taylor's series gives 


ft! 


0 


r = r o + 

o 

o 

+ (- M r Q / 

r + o ~r ) ip 2 /2 

0 0 0 

1 

4 

Of r r 

0 0 

ip 3 / 3! + a (_ 

* 

M "r _ / r + °\ r” ) 

0 0 0 0 

* 4 /4! 

4 

2 i* 

“ r o r o 

5._ , 2, 

ip / 5 ! 4 a (- 

- u r / r + <7 r ) 

^ 0 0 o o' 

0 6 /6! 


+ . 

4 

or, collecting coefficients of r and r , 
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t =[1 - M —( 0 2 /2! + a 0 4 /4! + a 2 ^ b /6! + a 3 0 8 /8 !+...)! r 
1 r 

0 

r ( 0 + a 0 3 /3 ! + a 2 0 3 / 5! 4 a 0 /7! + ...) 

0 

^ 4 t 3 8 

+ cr ( 0/2! + a 0/4! + a 0 / 6 ! + « 0 / 8 ! + ...) 

which gives the value of "? corresponding to any value of W . 

The results above are conveniently summarized by defining the 
transcendental functions 

s = ip / 2 ! + a ip / 4! + oi ip /'b l + .... 

0 

s = ip + a ip^ / 3! + Qf^ ^ ^ / 5 ! + <2 ^ / 7! + • • • * 

7 4 Z 6> 3 3 

s = ip / 2 ! + o 0 / 4! + o ip / 6 ! + c* ip / 8 ! + . . . . 

s = 0 3 / 31 + a 0 5 / 5 ! + a 2 0 ? /7! + a 3 0 79! + .... 

3 

which have the derivatives 


0 


0 


s 1 = o s, 

0 1 


i - 


r _ 


0 


1 


i _ 


3 2 

with respect to $ . In this notation 

' = V r 0 S l + %'z* <“i 

gives t for any value of 0 . Also 
r = t' = r s' + cr s' + p s' 


0 1 


0 2 


= rs+ (T s, + fis„ 
0 0 0 1 2 


and 


a = r' = r . s _' + a„ s,' + I* s,' = r A «s , + % s Q + Ms 


0 0 


0 1 


0 


Vo + ( ar o + " )S 1 
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gives r and cr corresponding to any value of 4>. 

It is also convenient to define the functions 
f = 1 - X s 2 /r 0 

g = r 0 S 1 + "O S 2 
which have the derivatives 
V -- - » . 2 '/ r 0 

= - X ‘S *o 

gl 1 r 0 S l' + Vg' 

= r + cr s . 

o c or 

However, the equation for t as a function of ip shows that 
g = • - t 0 - X s 3 

is an alternate expression for g. Differentiation of this equation 
shows that 

g' = t' - ns' 

= r - It S o 


is an alternate expression for g 1 . 

In terms of f and g, 

r = f r + g r 
0*0 

gives 7 corresponding to any value of ip. Differentiation of this 

equation with respect to ip gives 

r = f r + g r 
0 6 0 

and differentiation of "r with respect to t gives 

r = f r + g r 
0 B 0 


12 



where 


f = v = (- fi s 1 /r Q ) (1/r) 


= W> 


and 


g = g' <P = ( r - /i s ) (l/r) 
= 1 - fx s / r 


or 


t • 

g = g' tp = (r s_ + cr s,) (1/r) 
s 5 0 0 0 1 

= (r„s„ + cr. s ,) / r. 

0 0 0 1 

Another form of the equation relating t and ip above is obtained 
by defining 

y = r(r • r) - n 

which satisfies the identity 

• • 

y = r (r • 7 - 2 ^ / r) + p 


Then 


= a r + n . 


* ■ V r 0 S l + °0 S 2 + " S 3 


so that 


= ^ + “si + a „ s -> + M s„ 

0 0 3 0 2 3 

= */* + o' s + ( a r_ + ) s_ 

0 0 02 0 3 


t = t + r */> + a s + v A s . 
0 0 * 0 2 r 0 3 


Also, differentiation of this equation gives alternate expressions 

for r and o . 

/ 

t = r = r + a s, + v 

0 0 1 r 0 2 

t ,t = r , = a = cr s + *v s,. 

A A » A 1 


o o 


0 1 
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In addition 


<?' = (7 . r)' 

= y 

so that 


7 / * 7+7 *7 - r 7 


? + ■?•(- /i rVr ) 



- M 






as 


r 
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1.4 GENERAL. SOLUTION FOR COORDINATES 

m 

The formulas above may be used to calculate r, r from given 


values of r , 

7 and n , t Q} t. First 


r o = 

Vr • 

0 0 


a = 
0 

• 

r o* r o 


a = 

V'o - 2 " /r o 


are determined. Then the parameter ip and its transcendental functions 

S = 1 + a 4> 2 121 + o' 0 4 /4! + a 0/6! + . . . . 

0 

S 1 

s 2 = 

3 2 5 
ip + a ip /3! + a ip / 5! 

! + a ipf / 7 ! + .... 

2 4 2 
ip / 2 ! + a ip /4\ + a 

0 6 /6! + a 3 0 8 /8! + . . . 

S 3 = 

3 5/2 
^/3!+ Qf 5 ! + a 

0^/7! + a 3 0/9 ! + .... 

are obtained 

by solving the equation 


t = t 

+ r s + <t s 0 + MS- 

001 0 2 3 



for ip • Then 


r 


Vo + "o 5 1 4 


\x s 


2 


and 


£ = 1 - M s 2 /r 0 e - “ * * 0 > - " s 3 

£ = - C,/(r V 8=1- »V' 

give the final solution for the coordinates. 


r = 1 r o + 8 r o 

i ‘ ? 0 + 8 *0 


These are the equations for coordinates given previously by the 
author (1965). 
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The great advantage of these formulas is their complete 
generality. They are a modification of equations which were 
derived by Karl Stumpff (1947). Stumpff's derivation differs 
somewhat from that above and he normalizes his equations for 
computational convenience and simplicity. Similar formulas 
have been obtained independently by Herrick (I960) and Sperling 
(1961), but Stumpff was apparently the first to recognize that a 
generalized solution for coordinates of the two-body problem 
could be obtained by utilizing new transcendental functions sim¬ 


ilar to Sq, s j , s^, s above. Stumpff's work (1947) is given in 
his textbook (1959) and is also available in English (1962). The 
formulas above differ primarily in that the equations are not 


normalized but are left in a form that is continuous through 
the trival cases where p or (t - t ) is zero. The formulas are 
thus equally valid for negative values of p or (t - t ). 

The transformation from the variable t to the variable ip 
is often treated as a regularizing transformation since the sol¬ 
ution for r above is continuous through a collision which can 


occur if fi > 0 and 7 is parallel to 7 . For such a collision 

* u o 

r is indeterminate at the origin but is zero at the origin and 

continuous through the collision. Also, r and cr equal zero at 
the origin and are continuous through the collision* The con¬ 
tinuation of the solution through a collision is of course of 


little physical importance. However, it is practically advantageous 
because numerical problems are eliminated in the computation of 
position coordinates for near-collisions. 

The regularizing transformation from t to ip was used by 
Sundman (1912) in an investigation of the three-body problem. 
Stumpff realized that the transformation was of computational 
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value in solving the two-body problem, and that the resulting 
formulas could be placed in a form that was equally valid for 
elliptic, parabolic, and hyperbolic orbits including the cases 
of circular and rectilinear motion. He also realized the re¬ 
lation of the general formulas to classic formulas for elliptic, 
parabolic, and hyperbolic motion. This provides an alternate 
way of deriving the general solution, i. e. it may be obtained 
by manipulating classical formulas. Herrick (I960) used this 
approach and his formulas were modified to obtain the general 
solution above before Sconzo informed the author of Stumpff's 
prior work. 
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1,5 CORRESPONDENCE TO CLASSICAL FORMULAS 
In the case where a is negative, let 

-Jt- = - a 

a 

define the length a and let 
AE = y /M /a ip 

define the angle AE. Then the functions s^, s^, may be 

expressed in terms of a^and trigonometric functions of AE. 

For example, 

s = i/; 3 / 3! + a i^ 5 /5! + a 2 /7! + a tf 9 /9!+.... 

3 

= ^ 3 /3! - ( n /a) ^ 5 /5! + ( M / a) 2 ^ ? /7! - ( M /a) 3 ^ 9 /9! +... 

=y / / a 3 3 / 3 ! — V' q /a il> / 5 ! + y/ q / a £ 1 7 ! — . . . ■ 

\/Wa 3 

= ( n /a. ip ) / 3 ! — (y/ |i/a ji) / 5! + (\/ M/a ^ ) / 7 ! - ...♦ _ 

3 

v/ v /a 

= A 3 E/3! - A 5 E/5! + a'e/7 ! - .... _ 

y/T/a 3 

= AE - ( AE - A 3 E/3! + A 5 E/5! - AE/7! - -) 

y/lTU 3 

= A E - sin A E 

v^TTr - 3 

with similar derivations for s , s^, s^. 
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The result is that 


s o = 


cos A E 


sin AE 


v/ M /a 

1 — cos AE 
( M /a) 


A E - sin AE 


M 3 

These formulas give the following alternate solution when a 
is negative (which can occur only when p is positive). First 

y ? o • ? o 


r o 


a = 
0 


r ■ r ^ 
0 0 


(p/a) = 2 M/r Q - r Q • r Q 
are determined. Then solution of the equation 


t = t + r 


sin A E 


0 0 v/TT 

gives A E. Then 


r = r^ cos A E + & 


+ T 


0 


1 - cos AE 
( M/a) 


+ M 


AE — sin A E 




sin A E 


° s/ M /a 


+ M 


1 — cos AE 

( M /a) 


and 


f = 1 - 


0 


f = - — 


a 

r r 


1 - cos AE 
fi / a 

sin A E 


g = t - t() - p 


AE - sin AE 


i 

g 


= i - 


A 7r 

M 1 cos A E 


o v/m/ 


( M /a) 


give 


7 = f 7 o + 8 r o 


r = 1 r o + 8 r o 
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This solution is also valid when or is positive although AE and 
\/T /a are then imaginary quantities. However, the solution is 
not valid when a is zero and its numerical accuracy decreases 
as a approaches zero. 

The formulas for t, r, f, g, g above all reduce to classical 
forms such as those given by Herget (1948) if 

AE = E - E Q 

where E and E are, respectively, the eccentric anomalies at 
time t and at time t^. The well-known relations 

r Q = a ( 1 - ecos E ) 




a e sin E 


0 


are used to obtain the results.For example, 

sin AE 1 - cos AE 


t = t _ 4 r 


0 0 SU 


4 o 


0 


4 p 


AE - sin AE 


u /a 


'/Tj 


= t + a (1 - e cos E ) sin ^E / \/ /it /a 


J bC\x a e sin E (1 - cos AE) / ( n /a) 


4 fx ( A E - sin A E) / /7/a 
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so that 


s/v- 


/a (t - tg) = (1 - e cos E^) sin AE 

+ e sin E (1 — cos AE) 


+ AE - sin AE 


sin AE - e cos E^ sin A E + e sin E 

0 0 

- e sin Eq cos Ae 


+ AE — sin AE 

= AE + e sin E - e (sin E. cosAE + cos E 

U 0 i 

= AE + e sin E^ - e sin (E^ + AE) 

= (E - E ) + e sin E - e sin [ E^ + (E - E A ) 
0 0 D O 

reduces to the classic form of Kepler’s equation 


E - e sin E = E^ - e sin E + 

0 0 




/a 


(t t Q ). 


Since all the steps above are reversible, an alternate derivation is 
to define 

. E - E 
0 = 0 


\Z~li / a 

and express Kepler’s equation and the classic relations 
t* = a (1 - e cos E) 


f = 1 - * 

r 

0 


g = t - i 0 - 


[ 1 - cos (E - E ) | 
(E - E ) - sin (E 


V 




f = - 




r o 


sin (E - E ) 

0 


I 

g 


1 ~ — 
r 


f 1 - cos (E - E q )] 


sin AE) 


] 


21 



in terms of ip and ol . 

In the case where a is positive, let 

AF = \Z~a~ 4> 

define the angle A F. An equivalent definition 

A F = i AE 
where i = ri and 

a = - M /a . 

Either definition gives 
s 


0 


= cosh A E 
sinh F 


s i = 


a 


cosh AF — 1 
Of 

sinh AF - AF 


S 3 " 


j~d 


"wl^ich may be reduced as above to obtain classic formulas for 
hyperbolic orbits. The solutions in terms of A E and AF are 
both equally valid for all cases where a is not zero, but im¬ 
aginary quantities are eliminated by utilizing AE for negative 
a and A F for positive qj . However, both solutions decrease 

in accuracy as Ot approaches zero. 

In the case where a is zero, the functions s^, s^, s^, 


reduce to 
s 


0 


= 1 


s i = * 
s 2 = * 2 /2 

s = !/) 3 /3 ! 
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since only the leading term in these series is then non-zero. The 
resulting formulas may also be reduced to classical formulas for 
parabolic orbits. This suggests the possibility of using AE for 
negative a , ip for zero a , and A F for positive a , which 
corresponds to the classical approach. However, use of ip in the 
general solution for all values of Of is advantageous for several 
reasons. Computer storage requirements are reduced since only 
one formulation is required. The general soltuion is also contin¬ 
uous through zero & so that small values of a do no require a 
separate formulation for nearly-parabolic orbits. These advantages 
can become very important in practical applications where diff¬ 
erent types of orbits may be encountered. 
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2. THE PARTIAL DERIVATIVES OF THE SOLUTION 


2. 1 DIFFERENTIATION OF THE SOLUTION FOR COORDINATES 

The solutions r and r" are functions of the initial conditions ~r 
u 0 

and as well as /z and the times t and t . The differential re- 

» # 

lationships giving dr* and dr in terms of dr^ and d"? as well as d/i 
and the time differentials dt and dt^ are very important. The 
differential relation between d? and dr^, d - ?^ is given by Herget in 
a form due to Bower (1932) which is similar to the earlier work 

of Kuhnert (1879). Sconzo (1963) has extended this approach to give 

* . 

the differential relation between dr and dr^, dr as well. The 
approach given here is similar but concise expressions are ob¬ 
tained which are valid for all cases of two-body motion. The 
differential relationships are obtained by first differentiating the 
equations of the solution for coordinates and then combining the re¬ 
sults to eliminate all differentials other than dr", dr*' and dr" , dr" , 

0 0 

dfi , dt, dt . 

The differentiations are as follows. 

Differentiation of 


"r = f "r + g r 

0 B 0 

r = fr+gr 
0 6 0 


gives 


d? = f d^O + g df 0 
* ? o d f + 7 o d g 

d? = * d * Q + g d7 0 

+ d f + ? o d i. 
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Then 


gives 


and 


f=l— ns /r 
^ Z 0 


df = -d„ = 2 /r 0 - „ d s 2 /r 0 + » * 2 d r 0 /t o 

= [ — (f — 1) d r Q - ^ d - s 2 d M1 / r 0 


6 ' * - *0 - " s 3 


gives 


f° r g gives 


0 


- d t — du s„ 
0 3 

■' d! 3 


\i d 

s 3 - s 3 d M 

+ dt - dt^. 


; expression 



3 i + 

a s 

0 2 



d s i 

+ s i d r o + 

d 0 d s 2 + s 2 

da o 

d r 

0 

+ r o ds i + 

"0 d 5 2 + S 2 

d C7 0 


Similarly, 

i = - H Sj^r rj 

gives 

df = - d„ 3,/^rJ- „ d 5| /(rg + „ 3, 

= - 1 dr o /r o ' ‘ dt/r d S i^ r r o)‘ d " s X r r o) 

and 

g = i ~v s J* 

gives 

dg = - d^s^/r - n d s^r + dr/r 

= |-(id s 2 - (g- 1) dr - s d / r. 
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But the alternate expression 


a = (r s^ + ^ s )/r 

* x o o or 

gives 

d g = (r Q ds Q + s Q dr Q + a Q dSj + Sj d cr Q )/r 

- (r o s o + CT 0 s i )dr/r2 

= l s 0 dr 0 * 8 dr + r 0 dS 0 + a O d3 l + S l % 1 /r - 

The differentiation of the equation 

* = l 0 + r 0 S 1 + % S 2 + " S 3 

is equivalent to equating the two different expressions above for dg. 
This gives 

- n ds. - s d n + dt - dt = s. dr + r ds + ° ds + s 


1 0 0 1 


0 2 


r 0 dS l + "o dS 2 + “ ds 3 = - S 1 dr 0 - S 2 d % * s 3 dl ‘ 


+ dt - dt 


Similarly, differentiation of 


r = r 0 S 0 + d 0 s l + * s 2 


<- _ 

is equivalent to equating the two expressions above for dg. This gives 
[- n ds - (g-1) dr - s d p] /r 


= [.dr - g dr + r 


o ds o + a o ds i + s i d<T o ] 


dr = s Q dr o+ s x da Q+ s 2 d M 


+ r 0 ds 0 + C7 0 dS l + * dS 2. 
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The differentiation of 

s 0 = 1 + “ ^ 2 /2!+ a 2 0 4 /4! + a 3 0 8 / 6! + . . . . 

^ ^ ^ 

Sj = $ + a 0/3!+ a 0 /5!+ a 0 / 7! + . . .. 

s 2 = 0 2 / 2 ! + a 0 4 /4!+a 2 0 6 /6!+ a 3 0 8 /8! + . .. 

s 3 = 0 3 / 3 ! + a 0 5 / 5! + a 2 0 ? /7! + a 3 0 9 /9! + .... 

is conveniently expressed by utilizing 

2 ^ o / 

^/2!+ ^ of ip /4! + 3 Qf ip / 6! + . . . . 

0/31 + 2 a 0 / 5! + 3 a 2 0^/7!+.... 

0 4 /4! + 2 a 0/6! + 3 a 2 0 8 /8 ! + . . . . 

0 5 /5 ! + 2 a 0 7 /7 ! + 3 a 2 0 9 /9 ! + . . . . 


9 s 


0 


act 

as 


l 


act; 
8s _ 


3a 


8s 


9a 


Term by term differentiation of the series and combination of the 
results gives 


ds = 

0 

a s d 0 + 

"0 

8a 

d a 

ds = 

1 

S Q d + • 

9s i 

8a 

d a 

d l2 = 

8 J d 0 + - 

3 s 

2 

8a 

"da 

dS 3 = 

S d 0 +- 

/ S 3 

8a 

-da. 
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Finally, differentiation of 


r = (r • r ) 
0 ' 0 0 


1/2 


°b r Q * r 0 


V r o - 2 " /r 0 


gives 


dr 


1 


1 


0 


2 (r 0 ' r 0 ' 2 2r 0 dr 0 


= d "? • ■? / r 

0 0 0 


d cr = r • dr + r - dr 
0 0 0 0 0 


da = 2 r Q • dr Q - 2d fx /r Q + 2 dr Q /r o 
= 2 ( Mr # /r # 2 + d„ /r ). 
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2. 2 THE COMBINATION OF THE DIFFERENTIALS 

* 

The combination of the differentials above to give dr and dr 
as a function of drdr^, d /i , dt, and dt^ is simplified by using 
the relations. 

^0 = g ^ - g ? 

- -f r + f r 

and 

fg - fg = i 

which are derived in Section 2.4 below. 

Substituting the two expressions for r^ and r^ into the formulas 
above for d?" and dr* gives 

dr = f dr^ + g dr^ + *Jj d f + ? Q dg 

= f dr^ + g + (g ? - g ~r) df 

4 (-f ~r + f *r) dg 

so that 

dr = f d? 0 + g d7 Q + 7 (g d f - f d g) 

+ r* (- gdf+fdg) 

and, similarly, 

cT? = f d rr* + g d*?\ + r d f + ? d g 
0 6 0 0 0 s 

= f d? Q + g d? Q + (g r*- gr) d f 

+ (-f t + f "r) d g 
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so that 


The 

(gdf- 

Firstly, 

(g d f - 


dr = f dr" 0 + g dr^ + r (g d f - f d g) 

+ r^-g d f + f d g). 

I 

expressions (g d f - f d g), (-g d f + f d g) and 

f d g), { —g d f + f d g) are evaluated in order as follows. 


f d g) = g[-(f -1) dr Q - ^ ds 2 ~ s 2 d M 1 / r Q 

- i(s l dr 0 + r 0 dS l + V s 2 + s 2 d V 

= [-fSj -(f-1) g/r Q ] dr Q - f r Q dSj -(g H /r Q + f a Q ) ds„, 

- f s 2 dff 0 " S 2 /r 0 ) ^ 

9 s 

= [- f s i - t*" 1 * g^o 1 dr 0' fr 0 (s 0 d 0 d “ > 

0S 2 

-(g M /r Q + f Oq) ( 8l d^ + -55— da) 

-f s 2 da Q _ (g s 2 /r Q ) d M 

= [-fSj - (f - 1) g/r Q 1 dr Q - [ f r Q s 0 + (gn/r Q + f cr Q ) s ^ dip 

8S 1 - • 8s 2 

_lfr 0~9^ + (g /r Q + f cr Q ) -gjp 1 da 

- f s 2 d a 0 - (g s 2 /r Q ) d M 


= [-f Sj - (f - 1) g / r Q ] dr Q _ [f 


r S + O’ s MS 

0 0 0 1 ^ 


1 


r r 


g] r d ip 


MS 1 9s i 

- r- - r —-— 

1 r 0 da 


0 


+ ( 


r s„ + cr s 
0 0 0 1 M 


0 


MS 1 9s 2 

r r Q V 9a ^ d 01 


- f s ? d <r 0 - ( g s,/r n ) dg 


2 ' 0 
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[ ~ f S 1 ~ (f_1) g /r o ] dr 0 -<fg-fg)r 


d if) 


+ (s 

r 1 


® s i 3 s 

^ s 2 9 

”35 0 55~ ) d ot -fs^d cr ^ 


- (g s 2 /r Q ) 


= [-f 


*1 -(f- 1 ) g/r Q ] dr Q + _Jf_ 2 (s 


as 


as 


1 3 a; 


- s 


2 . do; 

5S - ) ~3 


~ f s 2 dCT o ■ ( 8 s /r ) dM 


But the identity 

2 


9 s 


S 2 = 2 < s ! 


1 


2 0 
3 s 


da ~ 5 

is treated in Section 3 . 1 . 

This gives 


0 ~da ) 


(g d f - f d g) = 


f-f s -(f- 1 ) g/ r ] dr 

x 0 i 

-Ms 


2 

r S 2 


2 ( M dr Jr + r • dir 

Oo 0 0 


d M /r Q )/2 


s , ( M dr /r + T . dr’ 
2 0 0 0 0 


- f s 2 d a _(g s^) d#l 

= f-f 8j -(f-1) g/r Q ] dr Q - (g_i) 

f S 2 d a 0 ' ( 8 S 2 /r 0 ) ^ 

= [-f Sj —(f— l)g/r Q +[(f-l) (g-ij/r^jd^ 

-(g-1) s 2 • dr Q - fs^ d a 0 - (s 2 /r Q ) dfi 

-- l-f »! -«-»/r 0 ] dr 0 - r Q /r 0 - (g- 1 ) ^ ? . d * 


d M /r ) 


fS 2 (? 0 * dr 0 + V d * 0 ) ~ {s 2 /r 0 )dfi 
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or, finally, 


(g d f - f d g) = - 


f s 1 + (f-l)/r 


0 


? • d?„ — f s r . dr 

0 0 2 0 0 


• • 


- ‘ S 2 V^O - ‘S- 11 S 2 V d? 0 - (S 2 /r 0> d " 


Secondly, 


(-g d f + f d g) = -g [-(f-l) dr Q - nds 2 - s 2 d/i ] / r Q 


+f fs l dr 0 + r 0 dS l + T 0 dS 2 +S 2 d V 


= Ifs + (f-l) g/r Q 1 dr Q + f r Q d Sj + (f a Q + g M/r Q ) ds 2 


+ f s 2 d o’ + (g s 2 /r Q ) dn 


= [f s, + (f-l) g/r ] dr Q + f r Q (s Q d ip + 


3 s 


da 


da ) 


d s 


+ (f + g M / r Q ) (Sj di ^ + 


— da) + f s_ d CT 
go; 2 0 


+(g s /r Q ) 


t f s i + (f ^ g//r o^ dr o + ^ r o s o + a o s i) + gM Sj/r 1 d ^ 


9 s ds^ 

+ [tr 0 “sr + (f "o* e " /r o )- 8Er ] d “ 


+ f s d<7 * (g s 2 /r 0 > dM 


[ f S + (f-l) g/ r 0 1 dr 0 + (f 


3 s 


r o s o + Vl 


9 s 


“ M s 


r r 


g) r d^ 


0 


+ [fr _i- 4 (f^ + gM/r ft ) -^] d« 


0 9a 


0 


0 9a 


+ fs 2 d o 0 + <gs 2 /r () )d„ 
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But 


= [f s 1 + (f-1) g/r Q ]dr Q + (f g - f g) rdf 


+ [f r 


9 si 9 

0 17 + (f °0 + 8 ' ,/r 0 ,_ Sr Jd “ +fS 2 d V (S W'*'* 


r as 

= [ fs 1 + (f-1) g/r 0 l dr 0 + rd*+| fr 0 — 


9 s 


-] 


+ (f a 0 + ^ /r o } da + f s 2 d a 0 + (g S -’ /r '' ) d(i 


2 0 


r O dS l + % dS 2 +,idS 3 = 


from above so that 


9 s 


r 0 (s o d4>+ “BeT da ] + 


- Sj dr Q ~ s 2 d % ~ s 3 dfl + dt ~ dt 0 


dS 2 0S 3 

( s , d * + TTZ .— a) +*t (s di^+ da ) 


0 ' 1 


9a 


9a 


= “ s l dr 0 - = 2 % - s 3 6 I‘ * dt - dt 0 


<r 0 S 0 + Vl =-< r 0 


3 s 3s_ 3s 

1 2 3 

- + a n - + M -) d a 

da u 8 a da 


s dr 
1 0 


S 2 d *0 " 


s d u 
3 H 


9 s 


9s 


+ dt — dt 


9 s 


0 


r dip = - (r 


0 9a 


+ a 


0 9a 


+ M—7—) da 
9 a 


s l dr 0 - S 2 d °0 - S 3 d “ +dt - dt 0 • 


This gives 

(-g d f + f d g) 


[ f s 1 + (f-1) g/ r 0 l dr 0 - ( r 0 


3 s 1 3s 3s 

- + <f. - +/i-) da 

da u da 8a 


- s dr - s d a - s_ dj* + dt - dt 


10 2 _ 0 '3 
3s 


0 


+ [f r 


1 


9 s 


0 


9 a 


+ (f a n + gfi /r Q ) -Ida 


9 a 


+ f S 2 d a 0 + < 8 S 2 ^ r 0 ) dfi 
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= [ (f-1) s + (f-1) g/ r 0 1 dr 


0 


8 s 


8s 


8s 


- r 


0 


- cr 


8a 


0 


8a 


- M 


+ f r 


8a 


0 8a 


9 S Z n 

+ (£ °o + *“ /r o> IT 


dad (£-1) d + (g s 2 / r Q “ 


0 


s ) djx 


-{- dt - dt 


0. 


But the coefficient of d a is 

8 s 

(f-1) r 


1 


0 


8a 


+ % + g “ /r o' to 


8 s 9s. 

2 3 

__ _ fX 


8a 


s 


1 


0 


= - M 


0 8a 
8 s 


+ t 


-jz s 


8 s 


0 


O' + (r s + CT s ) P/r 1- 

0 0 1 0 2 0 Ba 


(s 


1 


8a 


8 s 8 s 

- s -!) ♦ 3 

1 da 


8a J 


8 s 




8a 


and the identity 


S 2 S 3 = 2 (s 9S 1 s, 9s 2 ) 


- 1 


8a 


8a 


is treated in Section 3. 1. This gives the following coefficient for da 

9 s 9 S 3 

-j* (s s /2 + -)*= (g s, - g s o - /x - 1 X 2 — )/2 


8a 


2 3 
8 s 


8a 


= is =2 - s 2 (r 0 S l * Vz + 2 

a s 


3 2 


a« 


2 ’ ~2 ' 0‘ ace 


3 ]}/2 


Thus 


= (g s o - [ s , + M 2 


3 

(-gdf+fdg) = [ (£-1) s 2 +(f-l)g/r Q ]dr Q + U s 2 


8 s 


- [ s ? ^ t_t n^ + M 2 — 
2 0 ao- 


3 ) 


da / 2 + (f- 1) da Q 


+ (g s /r Q - s 2 ) dfx + dt - dt 


0. 
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Setting 


U 



gives 

(-g d f + f d g) = 


and, finally 
(-g d f + f d g) = 


as. 


(t-tg) + /x • 2 


da 


[ (f- 1) S 1 + (f- 1) g/r Q ] dr 


0 


+ (g s 2 - U) ( /x dr Q /r o + ? Q • df - d fi/ r Q ) 


+ (f-1) s 2 d o- Q + (g s 2 / r Q - s 3 ) dM + dt - dt Q 


[ (f-1) s + (f- l)g/r - n s 2 /r 0 )g /r 0 l dr 0 + U P /r 0 } dr 0 


+ (g s — U ) r • dr +{f—1) s d or 
6 2 0 0 2 0 


+ (-g s /r + U / r + g s / r 0~ S 3^ d M + at ~ dt 0 


[ (f-l)s 1 + (f-l)g/r Q - (f— l)g/ r 0 1 dr^- r^r 


0 


+ U (- g /r ) dr Q ‘ Tq/Tq + (g s 2 -U)r Q .dr o 


+ (f-1) s 2 d + ( U /r Q - s 3 )dfi + dt - dt 0 


[(f-1) s /r |?.dr to (-M 7 0 /r )• d7 0 -U ? 0 'dr 0 


+ g s„ ~r • dr" + (f- 1) s (r • dr" + "r • d7 ) 
® 2 0 0 2 0 0 0 0 


+ ( U / r — s ) d M + dt - dt 

' 0 3 0 


(f-1) s 


1 


r Q S Z r 0' dr 0 +S " 1)S 2 r 0' dr 0 


* • 


• • 


+ * s 2 r o' dr O + U r 0' dr 0 ■ u V dr 0 + < U /r 0 - s 3> d “ 


+ dt — dt . 
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The quantity U may be placed in the form 


where 


U = s (t — t )+M(^s — 3s) 

2 0 4 5 


^ 4 / 4! + a ^/(y\ + a ^ ^^/8!+ a ^ #*^/10!+.... 


S 5 = 


ip / 5 ! + a 
by using the identity 

8s i 

= * s 4 “ 


7/ 2 9, 

0 7 ! + a 0 /9! 


a 


0 11 / 11 ! + . . 


9a 


3s 


treated in Section 3.1. 

Thirdly, 

(g d f - f d g) = g {— f dr^/r^ — f dr/r - (j ds^ / r - dfi s^/ r r Q ) 

• , 

-f (s Q dr Q - g.dr + r Q ds Q 4 a dSj + Sj d a )/r 
= (-f g/ r Q - f s /r) dr - (f r Jr) ds -(gM/( r r o) + 


- (f s j / r ) d CT 0 “ (k s j/ r r Q ) d M 

f r 8s 

= (-f g/r - f s Jr) dr - - { a s d0 4 - dm) 

° 0 0 r 1 da 

9s i 

-(g H /(r r^)+ f a Q /r) (s Q d0 4 —— da ) 

-(f Sj/r) da Q - (g s l /(r r Q ) dfi 


= (-f g/r Q - f s Q /r) dr Q 4 [- f r^s^r 2 

f r 9s 

-(^/ryf yr)8 0 /r|rdM— 

8s i 

+ (g m/t r o)+ f Oq/t) —— ]d(* - (f Sj/r) d ‘ ^ S / f 
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But the coefficient of r is 


-f r 0 aSj/r -(g M /r r Q + f °o/ r ) s 0 / r = 


^ S i r o “ s i 


/r s + as, 
[ 0 0 0 1 


M s 


r r 


0 


r r 


0 


r r r 
0 

a \ 
— 

r ' 


0 


/ 2 

(s 0 ~ as i > 


and the identity 


2 2 

S — ot S 
0 1 


= 1 


is treated in Section 3. 1. 


Thus 


0 • 


(g d f - f d g) = 


(~f g/ r 0 “ f s Q /r) dr Q 


8s 


1 


8 s 


M— [ - (r - 

3 1 1 0 da 


a 


0 da 


ds 


+ M 


da 


) da 


- s l dr 0 - S 2 d % - S 3 d>1 + dt - "o' 


, fr o 8s o + ^ /r o tf ^0 


da 


ds 


da 


Bor 


-(f s 1 /r) d - ( g s^r r Q ) 


dM 


[- lj/>„ - f V + 




s l 1 dr o 


d s 


(r 


0 da 


3 s 


+ o’ 


0 da 


-f M 


8s 

Tor) 


fr 8 s 

0 0 


8 a 


g /i/ r 0 + f a 


0 


S 1 " 


8 a 


da 


37 



+ ( 




s - 
2 


f s 

— L ) d % + ( f 

r 0 3 


s _ - 


g s 


r r 


) dfx 


0 


(dt - dt 0 K 


But the coefficient of do? is 


A* 


3s 


(r 


0 8a 


+ O' 


0 8a 


r s + cr s 

I ( 0 0 0 1 


0S 


+ A* 


8a 




0 


M s 


■) + 


r r 


0 


/i s 


r r 


0 


r 9 s 

0 0 


0a 


0 s 


c ) 
0 


0a 




3 Ir o 


0 s 


+ (7 


3a 


0 s 


0 8a 


+ M 


0 S 


8a 


+ r (s 


0S 


0 


— s 


0 S 


8a 


0 


8a 


L ) i 


and the identity 

2 (s 


0 s 


0 


1 


- s 


1 8a 0 8a 

treated in Section 3. 1 below gives 


•> = S 1 S 2 * S 




r 


0 s 


1 


0 8a 
0 s 


+ cr 


9S 2 
0 0a 

0S 


+ M 


8S 3 S 1 S 2 +S 3 

- + r --- 


8a 


1 


3 [ 0 


_J_+ a Q U_+ M_!!3_ + r V2 + (r 0 s 0 + a Q s i + 

8a 8a 3a 2 ^ 


8 s 


S 1 S 2 

r -T^ + r n < - 

2 0 0Q/ 


S s ^ 0S 

L + + a 0 ( - 

2 0 da 


s s 
2 3 

+ H —-- + I* 


d s 3 "» 


2 'da 

for the coefficient of da . Also, the identities 


S 1 S 3 


9 s 


1 


+ s s = s s 
8a 0 3 12 


0s 


2 X 2 

- + S 1S “ s ? 

0 a 13 2 


treated in Section 3. 1 give 
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\ + -0 S 1 S Z + % s z + '■Vi *'*- 2 IT 1 72 


[s ( r s + cr s + g S„) + ^-2 -—3 /2 

£ U 1 U ^ J a™ 


- (g- 1) Sj/r + -^y [s 2 (t-t Q ) + g . 2 


] / 2 


= [ - ( g- 1)s i /r+ “3-U1/2 


for the coefficient of da • Thus 

(g d f - f d g) = [ -f g/r Q - f s Q /r + dr Q 


[ -(g- 1) Sj/r + 


M dr 


U ]( 


+ r • dr — 

0 0 


+ ( 


S 2 ~ 


fs gs 

-) d a- + ( -V— s - - ) d n 

r 0 3 3 r r 

r 0 


3 (dt " dt 0 ) 


= [-f g/r Q -f s 0 /r+ -t- s - (g-l)-!U 


U ] dr 


-(g-l)s i /r+ -y-Ul7 0 *dr o + ( 


— ) d a n 
r ' 0 


+ [ +(g- 1) 


r r 3 

0 r 


+ -tz S 

r n 3 3 

0 r 


] d M 


V (dt - dt o> 
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M-fg/r -£s /r + % s 4(g-l)f/r 0 + U |dr„ • r/:r 

r 


H— -JL 


+ ( - 


(g-D s 


1 


V ul r o' dr o + 


r r 


0 


, (g-n fs 
(- — dCT 


0 


o 


+ [ - 


M — / 
3 v 


r r 


0 


U 


S 3 ) ,dfl " ~3~ (dt - dt 0 > 


0 


M s 


[ - f s / r r + “ 

1 0 0 r r 


1 


0 r 


l 1 V d V ■V" | -'‘ ? o /r o »- d 'o 

r 

0 


(g- 1 ) s , . . f a +(g-l)/r ^ 

+ [- - + -4U 7 • dr"„ -- (i^-+ r^-drj 


0 0 


0 0 0 0 


+ [ - 


1 


- - ( 
3 v 

r r _ r 


0 


or ,finally * 
(gdf 


- f d g) = -f ( 


0 


r r 


0 


(g-1) s 


+ 


— + f* 


--s 3 ) ] d|i - -L (dt - dt Q ) 

r o 


1 1 


— > r 0 ,dT 0- -3~ U? 0' dT 0 

r o r 


• • 


3 U 1 r 0 dr 0" 


f s + (g-l)/r , 


(r „ • dr + r • • dr ) 
0 0 0 0 


+ [ - 


M , U 


r r 


- s^) ] d^t - Hj- (dt - dt Q ) . 


0 


0 


■ ♦ 


Lastly, since 

f g - f g = i 

it follows that 

d f*g + f d g - df*g - f dg = 0 

(- gdf + fdg) = -(g d f - f d g) 

f Sj + (f-l)/r Q 


» * 


0 


r • d r + f s ^ r*dr^ + fs^ r • d r 
00 200 20 0 


+ s z r 0 ‘d~r 0 + <(S 2 /r 0 )dli . 
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Then substitution of the final expressions for (g d f — f d g), 

(-g d f+ f d g), (g d f - f d g), (-g d f + f d g) into the expressions 

above for dr and dr gives the final expressions for the differentials. 
Firstly, 


d"r = f cf? q + g d "? + r* 


r f s + (f~ 1)/r 
1 1 ' 0 


0 


V “ ? 0 - ; =2 V d V 


* ^ 

2 ' 0' d '0 


- (£-1) s 2 r o -d? o - <s 2 /r 


o' •■<■] 


+ r 


(f-1) s 


1 


0 


V d *o + (£ -" s z V d ? o + <£ -" S 2 V d 7 o 


♦ # 


• # 


+ gs 2 V ? 0 + u 7 0' d7 0 - U V d "o + (U/r o' S 3 ) d ^ 


+dt — dt 

or, as the final differential dr, 


•] 


dr = f d r„ + U rr * dr 

0 0 0 


f s 1 + (f- 1 )/r 


0 


(f-1) 


r r -dr 

0 0 


£ =2 r r 0' dr 0 + r 


0 


S 1 r V dr 0 


* 4 


+ ( f -i) s 2 r r o .dr o+ g dr Q - U r r Q . d r 


• ♦ 


- £ s 2 r V dr 0 - 1b - £) s 2 ' V dr 0 + (£ -" s 2 r V dr 0 * g *2 r ~o' d *0 


+ [? (-s 2 /r Q ) + f ( U /r Q - s 3 ) ] dfi + 7 dt - ? dt Q 
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Secondly, 


• 9 


dr* = f d"? + g d? + ? 

0 0 


•01 1 . -► -► uit -► 

-f(- + — + —) r *dr - Ur • dr 

rr 2 2 0 0 3 0 0 

0 r r o 


H- 1 ' S 1 * A 




r * dr + —— U r * dr 
r 0 0 r„ 0 0 


f S + (g- 1 )/ r . 

1 . -H 


(r *dr + r *dr ) 

0 0 0 0 


°1 

-dpt - Hp ( 


r r 


0 


7 “ - s 3 } ^ (dt - d V 

0 r 


+ r 


f s 1 + (f- l)/r Q 


' 0' dt 0 + , *2 V dr 0 + ' S 2 V dr 0 


« » 


+ (g-1) T’o*dr" 0 + (s 2 /r Q ) dpt 


= f dT Q + u (- pt ~) ? Q -d? Q - f(^ 


0 1 1>\ 

+ — 0 + “o 1 r r « * dr,. 

r 2 2/00 

0 r r Q 


f s x + (g-l)/r f s i + (f-l)/r Q ^ ^ 


r r • d r + 
0 0 


0 


r r * dr 
0 0 


• * 


+ f s 2 ^? 0 -d? 0 + gd? 0 ) r Q .dT Q 

r 


(g- 1 ) s , . . f s + (g-l)/r 

1 -► “► 1 


r r *dr A 

0 0 


r r *dr 

0 0 


• # • 


+ f S r V dr o + (g-1) s 2 r r -dr +[r(- —) 

0 


+ f 


U 1 j* 

- f* -3 ( — - s 3 ) 1 M-5-( dt - dt ) 


0 
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or, 


w 

as the final differential cfr, 


dr = f dT + U 7 rl ' d"^ 


• s o 


ii f s +(g- 1)/r U 

-\~i + -2 + “ 2 > * ? 0 **0 - 7 - °* ° 

0 r r o 

fSj + tf-D/r ^ 

+ --- r r • dr + f s 2 r r *dr + g d r Q - U r V dr Q 

0 

_ f * d ^0 " (g - — 1 7 ^ 0 * ^0 + * S 2 77 0 d7 0 


• • 0 


+ (g-1) s 2 i r r' 0 *di r 0 +[ 'r(-s 1 /r r Q + 7 (s^r^ + 7 ( \J / r Q - s^) ] dfi 


00 00 

+ 7 dt - 7 dt 
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2. 3 THE PARTIAL DERIVATIVES 


The differential relations obtained above may be expressed in 

terms of a set of partial drivatives as follows. Let x } y, z be the 

# 

components of r and x, y, z be the components of "r so that x^, y , 
z Q are the components of r^ and x^, y Q , z Q are the components of 
Then the differential relations above may be expressed in the 
following matrix notation. 
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g 0 0 

+ 0 g 0 

0 0 g 



f Sj + (g-l)/r 
r 







The coefficient of dt in the above equations is a check on the solution 
of the differential equation. The coefficient of dt Q is the negative of that 
for dt which indicates that an increase in dt^ is equivalent to a decrease 

in dt and vice-versa. In most applications of the differential relations, 

• * • 

t and t are treated as fixed quantities and x, y, z, x, y, z are con¬ 
sidered as functions of x^, y , z^, x^, y^, z^ and sometimes ^ as well. 
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The differential relations are usually expressed in terms of a set 
of partial derivatives as follows. 



8x/9y Q 9 x/8z q dx 

dy/dy Q 8y,'8z Q dy Q 

dz/dy dz/dz dz 

0 0 _ L 0 J 


8 x/ 8 x q dx/ 8y Q 8 x/8z Q dx 

+ 8y/8^ 0 8y/8y 0 8y/8z Q dy Q 

dz/dx. dz/dy dz/dz dz 

u 0 0 J L 0- 




8 x/ 8y Q 
8y/8y Q 
8 z/ 8y Q 


8x/ 9z 
8y/ 8z 
dz/dz 



8 x/ 9 x q 8i/8y Q 8x/9z q "| j"di " 

+ 8y/8i Q 8y/9y o 8y/9z Q dy Q 

_8z/8x 0 9z/ay 0 8z/8z Q _ _dz 
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The partial derivatives of x, y, z and i, y, z with respect to y 
are merely the coefficients of d^ in the explicit differential expressions 

above. 


r 


r 


r ™ 

- S 2 /r 0 

- u/t o - V 

B x/ By I 
By/ Bp 


X X 

y y 

• 


B z/ By 

- 


z z 

_ — 




— 

9x/ By 


— 

• •• 

XXX 


- 5 i /tr V 

ay/ By 

- 

• * • 

y y y 


s 2 /r 0 

B^l By. 


* • * 

z z z 


U/ r - s 

L 0 3-1 


The partial derivatives of x, y, z and x, y, z with respect to 
x y , z Q and x Q , y Q , z Q are given by collecting the coefficients of 
dx , dy , dz and dx , dy A , dz A in the explicit differential expressions 

0*00 o 0 0 

above. A more compact matrix notation is also used to express the 
results as follows. 
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9 x/ 9x 

9x/9y Q 

3x/9z q 


f 

0 

0 



r n 

• 

X 

r n 

9y/9x 0 

9y/9y 0 

9y/9z Q 

= 

0 

f 

0 


u 

• 

y 

x o y o z o 

9 z / 9x 

9z/ 9y Q 

3z/ 3 z q 



0 

f_ 



« 

Z 

L. -J 

L _ 


y y 

z z 


r fs 2 + (£- 1) / r 0 


0 

(f-1) s 


1 


0 


— f s 


(f-l)s 


XV Z 
0 y 0 0 


x V z 
_ 0 y 0 0 . 


~3 x/ 3x^ 

9x/9y 0 

9x/ 9 z q 


g 

0 0 


X 

r 1 

9y/9^ 0 

9y/9y 0 

9y/9z 0 

= 

0 

g 0 

-u 

• 

y 

o 

• 

o 
• >> 

o 

•X 

9 z/ 9x 

9z/9y 

9 z/ 9z 


0 

0 g - 


• 

z 

L _ — » 


y 

z 


•n 

x 


y 

Z 


•-f s 2 -(g-l)s 2 
(f- 1)s 2 g s 2 J 


x o y o z o 


X 


0 


z ^ 

0 0 


3 x/ 3x Q 

9 x/9y Q 

9x/ 0z 


““ * 

f 

0 

• 

0 


« # 

X 

9 y/ 9 Xq 

9y/9y Q 

9 y/9 z q 

= 

0 

f 

0 

+ u 

y 

_9 z/ 9x^ 

9z/9y Q 

9z/9z q _ 


0 

L 

0 

• 

f_ 


«• 

z 


•• 

x 


0 


y 


o 


o 




x 

y 

z 


* 

y 


£ f- ♦ V + V > - 

r r o r r o 

f 8 + (f-l)/r 0 


0 


£ s 1 + (g- l)/r 


f s 


' x o y o z o 


L x o y o z o. 


3 x/ 3x^ 

9x/9y Q 

3x/ 3z 

9 y/ 9x q 

9y/9y 0 

3y/ 3z 

9 z / 9x 

L- 0 

9z/9^ 0 

3z/ z ( 


0 
0 

o-l 



g 

0 

o‘ 


mm 

X 


r 


0 

« 

g 

0 

— u 

• # 

y 


X* 

o 

o 


0 

0 

g 


«« 

z 


■* 



X 

v 


' f Sj + (g 

-l)/r (g- 



A 



r 






y 

• 

y 


• 



(g- 1) s 



• 


f S 2 





z 

z 







x o y o z o 
*o { o 


49 



In the expressions above for the partial derivatives, the parameter 
U = s 2 (t-t 0 ) + it (ip s - 3 s 5 ) 

is a monotonically increasing and unbounded function of the time t. In 
the case of elliptic motion where a is negative, all other quantities in 
the partial derivatives are periodic with the period 

T = 2 it \i jV -a 3 

of the elliptic orbit. In general, the partial derivatives of x, y, z , 

• * » 

y, z with respect to x^, y^, z^, x^, y^, z^ and ^ then have an un¬ 
bounded oscillation in time since the velocity components x, y, z and 
acceleration components 'x, y, z which are multiplied, by u have an 
oscillation about zero with the period T. However, this effect can be 
nullified for some derivatives if some of the other components x , y , 

Z 0 and V V *0 mult iplyi"g U are identically zero. In particular, the 
coordinate reference frame may be chosen so that only x and y are 

not identically zero, and the unbounded oscillation will then occur in 

only eight of the thirty-six partial derivatives of x, y, z , k , y, z with 
respect to x Q , y Q , z Q , i Q , y Q , z Q . 

In the case of elliptic motion, it is shown in Section 3. 1 that 


c - A E / 2! - (1-cos AE) 
4 ~ 


H /a 


s = A E/3! — ( A E-sin AE) 
5 ~ 5 ~ 


fyhich, with 


v/fT/T 


<P = 


AE 

M /a 


s = 1 - CQS ^ £ 

^ /a 

(t-O = r — m A E + (j 1-cos A E 


g = r 


\i /a 

sin A E 
J M/a 


/a 


, 1 - cos A E 

+ <T-— - 

0 

M /a 


A E - sin A E 
/—T“ 3 


M /a 
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1 


cos A E 


sin A E 


from above, gives 


U = 3pi 


AE - sin A E 

/ M / a 5 


p /a 




Of course, the partial derivatives may be computed for elliptic motion 
by using AE rather than ip to solve for the coordinates and determine 
all the quantities in the partial derivatives. Then the only secular term 
in all the derivatives is that due to A E in the expression above for U . 
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2. 4 THE INVERSE PARTIAL DERIVATIVES 

The solution of the two-body problem is still of course valid if 
t is treated as the reference time and r and r as initial conditions. 

^ m 

Then r and r are obtained as solutions of the differential equation 

0 0 

at the time t . For this inverse solution, 


a = r • r 

O' = r • r - 2 fx l r 

where OL is constant in time and therefore equal to the value above 

# 

computed from IT and . As in the solution above, If and its 
transcendental functions. 


= 1 + a f /Z! + a 2 ip 4 /41 4 a 3 ip 6 / 6 ! 


+ ... 


2x5 


= ip + a ip J / 3 ! + a c ip D fs ! 4 oc 3 ip 7 / 7 ! + ... 


s 2 “ ^ / Z \ 4 a ip / 4 \ 4 a ip /6!+a 2 $^/8!+.... 

s ="? 3 / 3 4 a ip 5 151 4 a 2 ^ 7 /7!+a 3 ^ 9 /9! + ... 
satisfy the equation 

l 0 = 1 + r ^l + al 2 + M V 

The bar is used to distinguish parameters of the inverse solution 
from those of the solution above. Also 


r 0 = r ‘o'" S 1 + “ s 2 


and 


f = 1 - /i s / r 

~ “ V r o r 


8 * ‘o - * -<“ 3 

1= 1 -K 1 2 /r 
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which give 


r = f r + a v 
0 5 


r = T~r + g r 

as the inverse solutions of the differential equations. 
However, 

1 r 

ip — / H t '{■ 

t 


f 

J r * 


o 


from the definition of ip above where r* is the value of the radius at 


the times t* between t and t. But 


ip = 


t J r^ 


dt* 


7 


dV* = - $ 


J* ''I s 


so that ip is merely the negative of ip . Thus 


s = 

0 0 


s i = - s i 


S 2 S 2 


s = — s 
3 3 


so that ip satisfies the equation 


t = t - rs, + as - 
0 12 


M s 


Also 


and 


r 0 = r S 0 - ""l + “ S 2 


f = 1 -M s /r 
= g 


g = t Q - t + M s 3 
= -g 


- (t-t - M s 3 ) 


f = M Sj/r r Q 


g = 1 — 11 s / r 

s 2 0 


= - f 


= f. 
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This gives the formulas 
r Q = g f - g t 

•V • - -V 

r = —f r + f r 
0 

which were used above to obtain the partial derivatives. In addition, 
substitution of 

*• 

r = fr+gr 

- i 4 .4* 

r = f r 0 + 8 r o 

into the first equation gives 
r Q = (f g - f g) r Q 

so that 

fg-f g = i 

which was also used to obtain the partial derivatives. Another property 
of the inverse solution is that 

U = s 2 (t Q -t) + M ( T 84 - 3 s*^) 

where 



which gives 

U = - U • 

The relations above may be used to express the partial derivatives 
in the inverse differential relations 
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dx ol f 8x 0 /8x 

dy o = 8 y o /8x 

UJ u 0 


^o^y a V 9z dx 3 x 0 /ai; 9x Q /ay 

%/3y 3y Q /3 z d y + 9y Q /aA 9y Q /ay 

8z 0 /9zj[_dzJ [_<V ai 3 z 0 /ay 




dip 9x / 9x 

dy 0 = a y 0 / 9x 

_ dz oJ l 9 V 9x 


9x Q /ay ax Q /9z 

3 y 0 / 3 y 8y Q /8z 

8 z /dy dz /9z 
0 0 


dx 

dy + 



ax Q /ax ax /ay 
3y 0 /9x 9y Q /9y 
az Q /ax az /ay 


9x /az 

3y 0 /9z 

az /az 




In particular, since 
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A similar method of expressing the partials of x^, y , z^ and 
x , y^, with respect to x, y, z and x, y, z by first using the 
functions of T and then converting to functions of 4 gives 


8 Xq / 3 x 

3x q / 3y 

9x / 3 z 


— 

0 

g 

0 

o n 


f— n 

X 

0 

8y Q /ax 

3y 0 /3 y 

3y Q /3 z 

- 

0 

* 

g 

0 

-U 

y 0 

3 z / 3x 

L 0 

3Zq/ 3y 

9z„/3z 

0 _j 


0 

w 

0 

* 

g^ 


>0- 


c 


] 


'0 


y 


o 


o 


0 "0 


aJ 


-fs^ (g-l)/r 
r 

_ (g-l) s l 


is z 

(g-l)s 


x y z 


• # 


x y z 


3x / 3x 

3x q / 3y 

3x / 3 z 


g 

0 

0 


x o 


9y Q / 

3y Q /3y 

3y Q /3z 

= — 

0 

g 

0 

+ u 

y o 

1—1 

3 z / 3x 

L. 0 

3z q / 9y 

SZq/ 3 z_ 


0 

0 



L*o- 






/— 

x o 


4 — 

x o 


~“ f S 2 

(f- 

U s 2 




— 

y o 

y o 


jd-l) 

S 2 

g s. 
2 





z 


m 

Z 









L o 


cu 









r— 







3 x / 3x 

3x q / 3y 

3x / 3 z 


f 

0 

0 


• 0 

x o 


3y o /3x 

3y Q /3y 

9y / 9 z 

y o 

= — 

0 

# 

f 

0 

-a 

1 • 

y o 

1-1 

x: 

3z / 3x 

— u 

3z q / 3y 

0 Z / 0 Z 


0 

0 

f 

f 


9 

« « 

t- Q J 




r s 

, o 

1 


1 

\ 

f S 1 + 

(f-D/r 0 n 


o 


x y z 

• • • 

x y z 


* • * • | 
y 


o 


o 


y 


o 


y 


o 


^ z o 




r r 


0 


0 


0 


f Sj + (g-l)/r 


f s 


x y 

• m 

x y 
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y 


9x^/ K 

3x 0 /9y 



f 

0 

0 


* • 

x o 

8y 0 /9i 

9y Q / 9y 

9y Q / 9z 

~ 

0 

f 

0 

+ u 

• • 

y o 

3 z Q /3x 

9z Q /ay 

9Zg/ dz 


0 

0 

f 


• • 

^ z 0- 



+ 



o 


o 


0 


f Sj + (f- l)/ r 0 


0 


— f s 


(f-l) Sj 


(£-1) s 2 J 



These inverse partial derivatives are often utilized for computation. 
However, they satisfy identically the well-known formulas 


3 x / 3x 

0 

9 X() /9y 

9x q /9z 

3 y Q /3x 

9y 0 /9y 

9y 0 /9z 

3 z / 3x 

0 

dz Q /dy 

9z / 9z 

3 Xq/ 3x 

9 x Q /9y 

9x / 9z 

9y o / 9x 

9y Q /9y 

9y Q /9z 

_9z / dk 

9Zq/ 9y 

9z /9z 

0 

9x / ax 

9x q / 9y 

9x / 9z 

9y Q / 9x 

9y Q /9y 

9y Q /3z 

9 z / 9x 

9z 0 /9y 

9z / 9z 

9x / 3x 

ax Q / 9y 

• 

9x q /9z 

9 y Q /9x 

3y 0 / 9 y 

9 ; o /9z 

9 Z^/ ax 

9z o /9y 

3z^/3z 


9i/9x 0 

3x/3y o 

3x/3i Q 

9y/9 x q 

9y/3y Q 

9y/3z Q 

Jtzj 9x q 

a^/9y Q 

3z/ 3i 

0 J 

3x/ ^ 

3 x/ dy Q 

3x/ 3z^ 

9y/ 3x q 

9y/9y Q 

9y/ 3 z q 

9z/ 9x 

9z/ 3y Q 

3z/ 3z 

0 J 

9x/ 9c 

9 x/9y Q 

3x/ az Q 

9y/ 3x q 

9y/3y Q 

9 y / 9 z q 

9z/ 8x 

3z/9y 0 

!>i/az 0 _ 

9 x/9x 

9x/ 9y Q 

— 

ax/ az Q 

9y/9x Q 

9 y / 9y 0 

9 y / & ? 0 

9z/ 9x„ 

L 0 

3z/9y Q 

3 z/ 3z 

0 J 


T 


T 


T 


T 


where the T indicates matrix transposition. Thus the inverse partials of 

Xq, y^, Zq and x^, y^, z^ with respect to x, y, z and x, y, z may be 

^ 0 0 

obtained directly from the partials of x, y, z and x, y, z with respect 

to , y^, z^ and x. , y , z . 

0 y 0 0 0 y 0 0 
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3.SOME USEFUL FORMULAS FOR DERIVATION AND COMPUTATION 


3. 1 THE PROOF OF SERIES IDENTITIES 

Several identities were used in the derivation above for the partial 
derivative s. 

As an example, consider 
2 2 i 

s o - “ 5 1 = ’• 

This may be proven by collecting coefficients of the same powers of $ 
in the products of the series. 

(l + ai l> Z /2\ + a 2 i/?/4! +_) • (1 + a ip Z /2l + a ip 4 /4l +_) 

-a (ip + a >P 3 / 3! + a Z j/ 5 /5 !+ ...)(if + o;!f 3 /3!+o; 2 if 5 /5!+_) 

O O /I 

= 1 + a ip (1/2 ! + 1/2 !) + a ^ (1/4 ! + l/(2 ! 2 !) + 1/4 !) + .... 

- a Z - a 2 i/- 4 (l/3 ! + 1/3 !) -_ 

= 1 + cn ip (1-1) + cl i/j 4 ( 1+6+1 — 4—4) /4! + .... 

= 1 


The product is valid since the series s . s , s , s , s , s above are 

0 1 2 3 4 5 

always absolutely convergent. However, an alternate proof is to first ob¬ 
serve that the identity is true for zero a where each series reduces 
to its leading term. If a is negative, then 

a = - {it a 


Sq = cos AE 
Sj = sin AE 

v/V /a 

From above where 

AE = / Jf a" 

so that 

2 2 2 u 

s^ — a s, = cos AE + —^ 

0 1 a 

= 1 . 


sin AE 2 
s/hTa 


2 

cos 


AE + sin 


AE 
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Similarly, if a is positive, then 

s = cosh AF 
0 


Sj = sinh’ AF 


\/a~ 


from above where 


AF ~\/a 4 


so that 


0 ^ 2 sinh. AF .2 2 2 

0 ~ a s \ = cosh AF -=—) = cosh AF - sinh AF 

s/oT 

= 1 . 

Thus the identity is true for all a and all . 

The other series identities used in the partial derivatives involve 


9s 9s 9s_ 9s ^ 

0 1 2 3 rni 

da 9 ~da 1 dcT 9 da * They may also be proven by collecting 

coefficients of the same powers of .ip in the products of the series. 

However, the relations 


0 = ip s 


* S 2 - S 3 


ip s - 2s 
3 4 


= ip s - 3s 
4 5 


may be proven by a term-by-term comparison of the two sides of each 
equation. If these formulas are substituted in the identities, each identity 
is valid for zero Q?. But 
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s 


0 


s 

s 

s 


1 

Z 

3 


1 + as 2 

ip + Of 

/ Z ! + Q! s 4 
^ ^ / 3 ! + as 

5 


similarly, so that, for non-zero a , 


s o- 1 


1 - S 


0 


s„ = 


a -a 

cos h AF — 1 1 — cos AE 


O' 




\i! a 
^ - s 

1 


s_ ~ 


a - Q? 

sin h AF/ \Z~cT - AF|V a AE/ v/ M /"a" - sin AE/ \/ pi /a 


Of 


^ /a 


sinh A F - AF 
3 

✓S' 2 

s - * / 2! 


AE - sin AE 

i/TTTa 3 

* /Z! - s_ 


a 


a 


_ (cos h AF -1)/ a - (AF/\^a) 2 /Z !_ ( AE/n/ 7T7^) 2 /2 !-(1 - cos AE)/( jx/a) 

*x/a 


a 


(coshiF -1) - AT/2! _ A 2 E/2! - (1 - cos AE) 

... - 2 


a 


( u /a) 


s - ^ / 3! 

s — J 

5 - 

a 


tp /3! - s 


a 


(sin h aF - AF)/ Vg •( AF// a~)^/3 !_ (A E^T/g) / 3-( A E- sin AE)/fy /o? 

a fx/a 

(sinh AF - AF) -A 3 F/3! AE/3! -( AE - sin AE) 


/aT 1 


/7J 
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Use of these relations proves the identities for both negative 
real or positive 0! where AF is real. For example 


a where AE is 


3s 


1 


3s, 


*2 s 3 = 2 ts 2~s; 

3s 


— s 


1 3a 


) 


1 

= s ( 2——) - s ( 2 
2 3a 1 


3s 


3a 


= s 2 ( ip s 2 - s 3 ) - s 1 ( ip s 3 - 2 s 4 ) 


is true when a is zero since 

(♦ Z /2!) ( >p 3 /3!) = ( <p 2 /2!)( t p //2! - ^ / 3 !)-( «HU # 3 /3! - Zip 4 /4!) 
ip / 12 = j/) 5 (1/4 - 1/12 - 1/3 + 1/12) 

= ip 1 12. 

For negative a , 


1 — cos AE 

\Zpj a 


AE - sin AE 



1-cos AE 
^ /a 


AiE 
W /a 


j.— c o s AE 
^ /a 


AE - sin AE ^ 
\/ /i / a ^ / 


sin AE j AE AE - sin AE 
v/V /a VM/a /m /a 3 

A E/2! - (1 -cos AE) \ 

u /*) 2 i 

(1-cos AE) ( AE - sin AE) = (1-cos AE) (AE - AE cos AE - AE + sin AE) 

2 

- AE sin AE( AE - sin AE) + 2 sin AE A E/2! 

-2 sin A E (1 - cos A E) 

= (1 - cos AE) (sin AE - AE cos AE) 

2 

+ AE (1 - cos AE) — 2 sin AE (1 - cos AE) 

= (1 - cos AE) (sin AE - AE cosAE + AE + AE cos AE 

- 2 sin &E) 

= (1 — cos AE) (AE - sin AE). 
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A similar result follows by using AF for positive a, but the use of 
A E is actually sufficient since the trigonometric formulas used are 
valid for imaginary arguments. 
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3. 2 THE USE OF DIMENSIONLESS SERIES 


For many purposes it is convenient to define the dimensionless 
argument 

A = a $ 

and its transcendental functions 

2 3 

c = 1 4 A/Z! 4 A / 4! 4 A /6! 4- 

® 2 3 

<5 ^ = 1/1! 4 A/314 A / 5 ! 4 A/714..., 

c^ = 1/2! 4 A/4! 4 A/6!4 A / 8! 4 . . . . 

c = 1/3! + A/5! + A/7! + A 3 /9 ! + . . . . 

c = 1/4! +A/6! + A 2 /8! + A 3 /10! + .... 

4 

c_ = 1/5! + A/7! + A 2 /9! + A 3 /11! + .... 

5 

which are related to s_, s,, s_, s , s , s by 

u I l j 4 o 



c, = 1/2! + A c 
2 4 

c, = 1/3! + Ac. 

D O 

or, conversely, if A t 0, that 
c & = ( c 3 ~ 1/3 !)/ A 

c . = (c, - 1/2!)/ A 
4 2 

c 3 = (Cj - 1)/A 
c 2 = < c 0 - 1 )/* 
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In terms of trigonometric or hyperbolic functions, for A ^ 0, 

2 3 

Cq = 1 +A/2! + A / 4 ! + A / 6! + .... 

2 4 6 

- 1 +/T /2! + /T / 4 ! +/T /6 ! +_ 

= cos h /*= cos h(i /- A) 

= cos y/ — A 

2 3 

c = 1 + A/ 3 ! + A / 5! + X / 7 ! + .... 

3 5 

= ( t/"x+ /T / 3 ! + /I / 5 ! + . ... ) 

= sin h /T //a" = sin h <! /Tx ) / (i /TT) = i sin \f-k ) 

= sin \/-X / /-A . 


Let Cq {4A} and c^{4A}be, respectively, the values of c^ and c^ for the 
argument (4A) rather than A . Then, forA>0 and A < 0, respectively, 
c q {4A} = cos h v 4A = cos 


= cos h {2\/ X ) = cos (2 /-X ) 

2 t 2 

= 2 cos h v A —1 = 2 cos v/ — A —1 
c^ {4A} = sin h V4 A / \/ 4A = sin V — 4A / l/ —4 A 

= sin h (2 /x ) / (Z\f\~ ) = sin (2 /- X ) / (2 \/- X) 

= 2 sin tt/T cos h/T / (2/x) = 2 sin ^- A cos v/ - X / (2 /TT) 

= cos h \/T~* (sin h.\f\~!\pk ) = cos y — A * (sin /T/ /Tx) 
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so that 


c 0 { 4 * 1 = 2 c o {x) ' 1 

c, < 4 X t= c 0 {Xt ' C 1 {X> - 

These two relations are also true for zero A , and are thus valid 
for all values of A . It is worth noting that the argument 

, ,2 

A = a ip . 

and a are simultaneously positive, zero, or negative except in 
the trivial case where ip an d thus A are both zero. 
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3.3 FORMULAS FOR THE SOLUTION OF KEPLER'S EQUATION 
It was shown above that the general equation 


relating t and 4> is a generalization of Kepler's equation for elliptic 
motion. Stumpff (1962) calls it the main equation but prefers an 
equivalent of the alternate form 


' = *0* '0 * * Cr 0 S 2 + * 0 V 

The first, second, and third derivatives of this alternate form 
with respect to & give the formulas 


r = 


r o + 


cr s + V s 
0 1 y 0 2 


for the radius, 

cr = cr s + y s 
0 0 *0 1 

for the scalar product of the position and velocity vectors, and 

y = r o s o + * °o S 1 

for the parameter y which equals n for parabolic motion, i. e. for 

01 ~ O' When CT 0 and y Q are both zero, r is always equal to r 

whereas a and y are both always zero. Then the motion is circular 
and f is a linear function of t. 

In computing coordinates with or without partial derivatives for 
a particular time t, it is necessary to solve for that value of 

which satisfies the main equation above. With a slight change in 
terminology, let 


(t - 


V = 


0 1 


+ CT 


0 2 


H s 
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be the time interval computed directly from any value of ip , 
and let t be the desired time interval. Then the change in 
ip required to change (t - t^) to t is given by 

T = (t - t ) + d * * V Alp 

° d^ 

to first order accuracy. Thus 
r — ( t - t ) = r 

where r is the radius for the particular value of ^ . The value 

ip* = ip + Alp 

is then a better approximation for the desired value of $ if the 
first order approximation is accurate. Explicitly, if 


At 



s i + 


or + 
0 2 



T 


1 = r oV Vl + “ S 2 


are first computed from then 

ip = ip - At/ r . 

This is Newton’s method for the solution of the main equation. 

Using *l> as a new value for ip then gives a new approximation, 

etc. However, a method of stopping the iterations is necessary. Also, 

Newton's method is not always convergent and should be backed up by 

* 

alternate methods for obtaining a better approximation ip from ip . 
The iterations may be controlled by always bounding the solution by 
the maximum known value ip and minimum known value ip + ior which 
the respective residuals At_ and At + are, respectively, negative and 
positive. Then any approximation ip is accepted at a new value of ip 
only if 

>P_ < < 4> + 

i. e. if 4> is between ip and ip + but not equal to ip_ or ip + . 

When the value $ computed by Newton's method fails to lie 
between and $+ , alternate methods may be tried to obtain a 
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satisfactory value of $ . One such formula is 
ip* = ip (1 - A t/t) 

which gives smaller corrections to $ for smaller values of I At 
A nother is to define ip as the value satisfying the interpolation 


formula . 

- 4 >* = 



At_ A t + 

Explicit solution of this equation gives 



At_ 

At + — At 



The use of these formulas is treated in the description of the 
Fortran subroutine in Section 4.1. 
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3.4 FORMULAS FOR SERIES COMPUTATION 


Each iteration in the solution of the general form of Kepler's 
equation requires the evaluation of s^, s^, s^, and from a given 

value of ip . The evaluation of s^ and from the final value of ip 

is also required for the determination of partial derivatives. However, 

it is convenient to compute the dimensionless series c^, c^, c^, c^, 

c , c of the argument 
4 5 

,2 

X = a ip 

rather than s rt , s., s_ , s . s . s directly. Then 

0 1 2 3 4 5 



so that s and are never required explicitly. 

4 5 

The series c , c . c , c . c . and c must be^computed from 

0 1 2 3 4 d 

arbitrarily large values of | aJ , i. e. for arbitrarily large values 
of |of|and |^| , in order to determine coordinates and partial de¬ 
rivatives for arbitrarily large values of |t|, i.e. for any value of t . 
This may be done by repeatedly dividing X by 4 when x\ > 1 until 
the m division reduces the argument to one or less in magnitude. 

Then c and c, are computed as described below with the reduced 
0 1 

argument. Then the formulas 
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c 0 = 2Cq {a} -1 

Cj (4 A} = c Q {A} * Cl {A} 


are applied repeatedly m times to obain their true values for the 
original value of A . Then the formulas 



C 3 = ^ C 1 ” 

c 4 = (c 2 - 1/2) / A 

Cg = (Cj. - 1 / 6) / A 

above give the true values of the four remaining * series. This approach 
seems to be quite satisfactory, but an alternate method is to determine 


Cq = cos V-A 

c 1 = sin V -A / V~“A 
directly when Ac.- 1 or 
c^ = cosh 

c ^ = sinh yfX / yT\ 

directly when A > 1 . Then the same four formulas above give c , c , c 

2 3 4 

and c ■ 

D 

The series c^, c^, c^, c^, c^, c^ may all be determined from a 
given value of A by first using the nesting formulas 



X 1 

* . . . !\ * 1 

(2n+3)(2n+2) 

_ 

(2n+ 1) (2n) 

) 9* 8j 

x 1 

* o 

(2n+2)(2n+1) 

(2n)(2n-1) 

J 8-7 



to determine c^ and c . The number of terms n required for the 
computation is considered below. Eight terms, i. e. seven leading 1 1 s, 
are sufficient for accurate floating-point computation to sixteen decimal 
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digits when | x\ < l which is always true since if | X\ > 1 initially 
it is either reduced as described above or the alternate method is 
utilized. Then the resulting values of c^ and c^ give 


c^ = 1/6 + X c^ 

c 2 = 1/2 + x c 4 

C 1 1 1 + Xc 3 
c 0 ' 1 * * c 2 


for 

and 


the 

c i 

c i* 


remaining four series. If the argument has been reduced, c 
are then used as described above to obtain the true values of 

c„, c„. c„, and c_ 

2 3 4 5. 
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3. 5 THE ACCURACY OF THE SERIES COMPUTATION 

The limiting cases in the nesting computation above occur when 
I a| = 1. If A = -1, then 

c = cos V -A = cos (1) 

= .54030 

C 2 = (c Q - 1) / A = (. 54030 -1)/ (- 1) 

= .45970 

C 4 = (c 2 - 1/2) / A = (.45970 - .50000) / (-1) 

= .04030 
and, similarly, 

c x = sin yT-\ /v/ — A = sin (1)/(1) 

= .84147 

c 3 = ( Cj -1) /A = (. 84147 - 1) / (-1) 

= .15853 

c 5 = (c 3 - 1/6) /A = (. 15853 - . 16666) / (-1) 

= .00813 

Also, if A =1, then 

c = cosh /A = cosh (1) 

= 1.54308 

C 2 = {c 0 " 1)/A = t 1 * 54308 - 1 )/ (!) 

= .54308 

c 4 = ( c 2 - l/2)/A = (.54308 - . 50000)/( 1) 

= .04308 
and } similarly, 

Cj = Sinh tTa / \Ta = sinh (1)/(1) 

= 1.17520 

C 3 = (c i " A = (1- 17520 - 1)/(1) 

= . 17520 

C 5 = (c 3 * 1/6)/ A = (. 17520 - . 16666)/(1) 

= .00754 
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Suppose that eight terms are chosen to compute c^ and c,. for 
I A.) = 1. Then if A = -1, the truncation error Ac^ introduced in 

computing c is 

4 8 9 10 

Ac = A /20! + A /22! + A / 24! + .... 

4 


= + 1/20! - 1/22! + 1/24! + - 

so that 

-18 

Ac < 1/20! = . 41103 x 10 

4 

from the properties of alternating series. 

Thus 



. 41103 x 10 
.04030 


-18 


10 . 


19 x 10 


-18 


Ac 


< .1019 x 10 


-16 


Similarly, the truncation error Ac,, introduced in computing c,. i 

8 9 10 

Ac — A / 2 1 ! + A / 2 3 ! + A /25 !+..•• 

5 


= 1/21! - 1/23! + 1/25! - - 


so that 


A c_ < 1/21! = . 01957 x 10 

5 


-18 


and 


Ac 


5 .01957 x 10 


-18 


. 00813 


= 2.407 x 10 


-18 


Ac 


5 1 2407 x 10‘ 17 
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If A = + lwith eight terms, the truncation error Ac introduced 

4 

in computing c^ is 

8 9 10 

A c 4 = A /20 ! + A /22! + A / 24! + . . . . 

= 1/20! + 1/22! + 1/24! + . ... 

< 1/20! + 1/(20 • 20!) + l/(20 4 • 20!)+ _ 

= (1 + 1 / 20 2 + 1 / 20 4 +_)/ 20 ! 

= [ 1 /( 1 - 1 / 20 )] / 20 ! = ( 20 / 19)/20 ! = ( 1 / 19!)/19 
= .04326 x 10" 17 

so that 


Ac 


.04326 x 10 
. 04308 


-17 


= 1. 004 x 10 


-17 


Ac 


. 1004 x 10 


-16 


Similarly, the truncation error Ac 5 introduced in computing c is 

8 9 10 

Ac 5 = A /21! + A /23! + X /25! + . ... 

= 1/21! + 1/23! + 1/25! + _ 

< 1 / 21 ! + 1 /( 21 2 . 21 !) + 1 /( 21 4 . 21 !) +_ 

= (1 + 1/21 2 + 1/21 4 + . .. . )/2l ! 

= [1/ (1-1/21)] / 21 ! = (21/20)/21! = (1 / 20 !)/20 
= .02055 x 10" 18 
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so that 


Ac 5 . 02055 x 10~ 
c 5 .00754 


= 2.72 x 10 


-18 



.272 x 10 


-17 


Thus eight terms are sufficient for sixteen significant figures 


in the nesting computation above for c and c^ with |A|< 1. A 
similar derivation will determine a sufficient number of terms if a 
different number of significant figures is required. Ordinarily this 


is determined by the particular electronic computer utilized, e. g. 
the IBM 7094 has 54 binary bits in a double-precision fraction which 
gives slightly more than 16 significant decimal digits for floating¬ 
point computation. 
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4. THE DOUBLE-PRECISION FORTRAN 4 IBM 7094 SUBROUTINE 


4. 1 DESCRIPTION OF PROGRAM INPUTS AND OUTPUTS 

The Fortran 4 IBM 7094 double-precision subroutine TWOBDY 
has the following nine inputs. The six initial position and velocity- 
components at the reference time t are input as a state vector 





* 



The seventh input is the time interval 


r = t - t 

0 

between the reference time t and the time t at which a solution is 
desired. The eighth input is the constant fi in the differential equation 
where 

M = G (m ^ + m^) 

if relative coordinates between two gravitationally attracting masses 
m j an d are being obtained. The ninth input is an approximate 
value for the final solution ip of Kepler's equation where 
♦ = 0 

is used if no approximation is available. 

The input value of ip is changed by the program to the actual 
solution of Kepler's equation for the input time interval r . Thus 
the ninth input to the program is also an output quantity. The number 
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of iterations required to solve the general form of Kepler's 
equation can be considerably reduced if a good input approximation 

for ip is available. In particular, solutions may be required for a 
whole series of values of t which differ by given increments. Then 
the use of the previous value of ip as an input approximation for ip 
at each new incremented value of T is very efficient. It is usually 
better to use an input approximation for ip even if it is only known 
very roughly. If not, it is more efficient to use an input value of 
zero although the program will work with any input value for the 
parameter ip . 

The program outputs the six position and velocity components 
in the form of a state vector 



for the time t. The 36 partial derivatives of x, y, z, x, y, z with 
respect to x Q , y Q , z^, x Q , y^, z Q are output in a 6 x 6 matrix 

9 x/ax Q 9x/9y Q 9x/9z Q ax/8i Q 9x/9y Q dx/dZ Q 

dy/ 8* 0 9y/9y Q 9y/9z Q 9 y/ 9 * 0 9 y/ 9 y Q 9 y/ 9z 0 

3z/ £bc 9z/9y^ 3z/ 3z^ 9z/ 9*^ 9z/9y^ 9 z/9z^ 

3x/9x 9x/9y Q 9x/ 9z Q 9x/9x Q 9x/ 9y Q 9x/9z Q 

9y/9x 0 9y/ 9 y^ 9 y/ 9z 0 9 y/ &Xq 0 y/ 

3z/3Xq 9z/9y^ 9z/9z^ 9z/9x^ 9z/9y 9z/9z^ 
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Similarly, the inverse partial derivatives are output in a 6x6 matrix 


9 x / 9x 

9x q / 9y 

9 y Q /9x 

9y Q /9y 

9z / 9x 

9z Q /9y 

9x^/9x 

9x Q /9y 

9y Q /9x 

dy^/dy 

0Z^/0X 

9z Q /9y 


9x / 9 z 

9x^/ 9x 

9y Q /9z 

9y Q /9i 

9z /9z 

9z / 9x 

9x^/9 z 

9x^/9x 

9y Q /9z 

9y Q /^ 

9z /9z 

9b / 9x 


9x Q /9y 

9 x Q / 9z 

9y Q /9y 

9y Q /9z 

9z q / 9y 

9z / 9z 

8* 0 / a Y 

9x^/9 z 

9y Q /9y 

9y / 9z 

y o 

9z Q /9y 

9z / 9z 
0 


The two 6x6 matrices above are also matrix inverses. That is, their 


product should give a 6x6 identity matrix. 

The six partial derivatives of x, y, z, x, y, z 
M are output as a vector 

9 x/ 9p 
9 y/ 9p 
9 z/ 9p 
9 x/ 9p 
9 y/ 9p 
Jz/ 9p_ 

Similarly, the six inverse partial derivatives of x , 
with respect to P are output as a vector 

9x q / dp"' 

9 y Q /dp 

9 z / 9p 
9 x /9p 
9y Q /9p 

_9 z Q /9Pj 


with respect to 


< • • « 
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Additional outputs are the acceleration vector 



at time t and the acceleration vector 



at time t A , as well as the radius 

0 tt . —2—r 

r = Vx+y+z 


at time t and the radius 



at time t^. 

The Fortran symbols for the inputs and outputs above are 
listed in comment statements at the beginning of the subroutine. 

The program may be used as a "black box 1 * by those not interested 
in the formulation or computation. However, a knowledge of the 
operation of the subroutine is very important if modifications are 
desired for particular applications. For example, it may be de¬ 
sired to eliminate the partial derivatives or at least make their 
computation optional since they are not required in many cases. 
However, such modifications will be worthwhile only where the 
particular application is encountered repeatedly. Ordinarily it is 
better to use the program as it stands without modification although 
this may involve some redundant computation. In any case, the 
operation of the program is described below for those who are in¬ 
terested. 
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4. 2 THE INITIAL COMPUTATIONS 


The following initial computations are performed only once 
before starting the series computation below that is repeated for 
each iteration in the solution of Kepler's equation. The following 
technique for computing the radius r Q at time t is used to avoid 
utilizing a square root subroutine. First 

b = maximum ( x 0 |, | y Q |, |z |) 

c = <V b)2 + ( V b)2 + < z 0 /b > 2 

are computed and the initial value of d is set equal to two. Then a 
better approximation d* for the square root of c is obtained from 


d* = (d + c/d)/2 

and d* is used as a new value of d to iterate the formula. The 
iterations are continued until d is no longer greater than d* which 
is then accepted as the square root of c. Then 

r ~ b • d ^ 

o 

gives the radius r^ at the reference time. The quantities 


cr 

0 


+ V Y + z 
0 0 y 0 y 0 0 





are then computed. It is also necessary to initialized a counter m 
to zero in order to compute the series described below for large 
arguments. 


In the trivial case where t is zero, ip is set equal to zero and 
the series computation below is performed. The known value zero for 


ip satisfies Kepler's equation identically so that all the outputs will 
then be computed by the program. If r is negative, it is known that 
ip must also be negative. Thus the right bound 4 : for >p is set equal 

T 

to zero and its residual At is set equal 
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to the known value - T , whereas the left bound 4>_ and its residual 
At are both set to i.e. the negative of the largest number 

available for computation. If t is positive, it is known that ip must 
be positive. Thus the left bound ip _ for ^ is set equal to zero and its 
residual At is set equal to the known value -r , whereas the right 
bound ip and its residual At are both set to + », i.e. the largest 

number available for computation. 

If the input approximation for ip then lies between *P _ and ip 
it is used in the series computation below. If not, Newton’s method 
using an initial ip of zero gives the approximation 

*= T /r 0 

which is then tried. If this value is not numerically between ip _ and 

* V 

ip - r 

is used to start the series computation. Thus an initial approximation 
for ip is always obtained which is between ip _ and $ + but not ec l ual 
to either. 

4. 3 THE SERIES COMPUTATION 

The computation of the series is begun by computing the argument 

X = a y 

of the simensionless series. If | x| > li the value is saved as ^ 
but X itself is repeatedly divided by 4 until |M 5 1, and the number 
m of required divisions: is also counted. Then 3c & and c 4 are computed 


by the nesting formulas 
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which in turn give 


c 3 = [1/2 + x (3c 5 ) ] / 3 

c. = 1/2 + Xc 

i. 4 


Cj = 1 + Ac 


c 0 = 1 + A c 2 

where each series is correct to sixteen significant figures. 

If the original value of X has not been reduced, the series s 

3’ 

s 2’ S 1 are com Puted as described below. If the original X has been 
reduced, i.e. if m> 0, the formulas 

* 2 

c o - 2 c o - 1 
❖ 

C 1 = C 1 c o 

jjc 

are applied m times with c and c as new values for c and c . 

5jC U ft t 0 1 

The final values of c and c are accepted as the frue values of c 

U 1 Q 

and c for the original value of A, i, e. for \ . Also 

P 

C 2 = (c 0 - 1)1 \ 
c = (c. - 1)/ \ 

J l P 

C 4 = (C 2 ' 1/2)7 A 

P 

3c = (3c, - 1/2)/ X 

D 5 p 

are computed as the true value of the dimensionless series. 

The dimensionless series c , c , c, give the series s , s , s 

5 4 1 3 2 1 

by the formulas w 


S 3 = 


4? c 


S = Ip 
2 


s , = 4> c 
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which are used to compute the residual in Kepler's equation below. 
The value of c is equal to s Q and will be used in Newton's method 
for solving Kepler's equation. If this particular series computation 
should be the last of the iterations for solving Kepler's equation, c 4 
and 3c will be used below to compute partial derivatives. The series 
computation described here seems to be quite satisfactory for ob¬ 
taining solutions for two-body motion. The absolute error in the so¬ 
lutions will increase as r and thus A increase, but this is an 
inherent error since the absolute error in ^ itself increases as it 
becomes larger if it is represented as a floating point number. 


4.4 THE SOLUTION OF KEPLER'S EQUATION 

The residual At for solving Kepler's equation is found by first 
computing 

8 = r 0 S 1 + g 0 S 2 
for the current value of 4> . Then 

AT = g + (1 S 3 - T 

gives the residual. The radius 

r = r 0 C 0 + Vl + MS 2 

for the current value of ip is also determined. If At is zero, the 

coordinates are then computed as described below since Kepler's 
equation has been solved. If At is negative, it replaces the old 
value of At and ip replaces the left bound ip_. If At is positive, 
it replaces the old value of At + and ip replaces the right bound ip + . 

Then Newton 1 s method 

$ 

Ip = ip - A T /r 

is used in an attempt to obtain a better approximation for tp . 

If ip* > ip and ip* < 4> , 4> is accepted as a new value of <P 

and the series computation above is repeated. Otherwise, the following 
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methods are used succesively in an attempt to obtain an approximation 

between ip and ip . First 
- *f 

4> = ip_ (1 - 4 At / r ) 


is tried if |Ar - |<| At+| but 
ip = ip (l - 4A r /r) 

“T *|> 

is tried if |Ar-| > At + I 


The factor 4 helps to give significant 

corrections to ip, in floating-point computation when |at - | or |at + I 

* 

or ip fails to lie between ip and ip , 

+ 


is small. If |At~| = |At + 
the value 


= 2 ip_ 

is tried if t > 0 , but 

ip* = 2ip + 

is tried if r<0 . This insures that the solution will always be isolated 
between finite bounds. If this method fails, the interpolation formula 
ip = - ( ip + - ip_) At/ (At — At ) 

I® tried. If this fails, the rnid-point 

>P = 4>_ + (4> + - >p_)/2 

between 4>_ and ^ + is tried. If this fails, the iterations are terminated 
and the coordinates are computed since a better value for $ is not 
obtainable. 


4. 5 THE COMPUTATION OF COORDINATES AND ACCELERATIONS 

The final values of the radius r and the function g have been 
obtained in the last iteration of Kepler's equation above. The functions 
Sj and s^ from the last iteration are used to determine 

(f-D = - ps 2 /r Q 

f = _ p Si / (r r o ) 

(g -1) = - H s 2 /r . 
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Then the coordinates are obtained from the relations 


X = x Q + [(f- 1 ) X Q + g X Q ] 

y = y Q + [ (f-1) y Q + g y Q ] 

z = z Q + [ (f - 1) z Q + g Z Q ] 

X = [f X Q + (g-1) X Q ] + X Q 

y = y 0 + (g- 1 ) y Q l + y Q 

z = [f z 4 (g-1) z ] + z 
1 0 6 0 0 

which give maximum accuracy with small values of t for which (f-1), 
• ^ 

g, f, (g-1) are each close to zero. Then 

/ 3 

x = - jxx/ r 

y = - fiy/r 3 

3 

z = - jiz/r 

give the acceleration components at time t, and the acceleration 

components 

.. . 3 

x„ = - y x / r 
0 ^00 

y o = •" y o /r o 3 

z o - - 11 z o /r o 3 


at time t are also computed. 

4. 6 THE COMPUTATION OF PARTIAL DERIVATIVES 


The partial derivatives are now computed by determining 

U = s 2 (t-t 0 ) + M (c 4 - 3c & ) Ip 5 
and using locations in the matrix of partial derivatives to store 
elements of the matrices on the right side of the formulas 
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— 

3 x/ 3 p 


r 

X X 


— — 

- S 2 /r 0 

9 y/ 3 p 


y y 


U / r - s 

3 z/dp 


z z 

La -J 


L 0 3 J 


F— *" 1 

3 *-! dp 


l '*~ • • • ~' 

xxx 


l 

1 

Cfi 

o 

9 y/ dp 

= 

» 

...y' y y 


s_ / r 

3 z/3*i 


• # # 

z z z 


o 

r\j 

VJ. 




L u/r o- s 3 J 


— —i 

3 x.~/dp 


X X 

0 0 


’ S 2 /r 

3 y Q /9^ 

— 

O 

O 


-U/r + s_ 

3 z /dp 


z„ z„ 


L. 3 — f 

L- 0 -J 


L o o J 



3 x Q /3/u 
3 y Q /3/i 

- 

n . - n 

x o x o x o 
y o y o y o 


r— "" 

V (r r 0 > 

s./r 

_9 Zq/9^ 


• • # 

Z„ z„ Z 

L o o oj 


2 

-U / r + s 

L 3 -J 


which give the partial derivatives with respect to ^ . Similarly, 
locations in the matrix of partial derivatives are used to store 
elements of the first matrix on the right side of the formulas for 
the 2x3 matrices 
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which are stored in locations of the inverse partial derivative matrix 
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Then the computations 


"3 x/ 9x 0 

9x/ 9y Q 

a y/ 9 x Q 

9y / 3y 0 

_9 z/ 9 x q 

9z/ 9y o 

~9x/ 9x 

9x/8y Q 

9 y/9 x q 

9y/ 3y Q 

_3 z/ 9x^ 

9z/ 9y Q 

3 x/ 3x^ 

9*/ 3y Q 

3y/3x Q 

9y/9y Q 

3 z/3 Xq 

9z/ 9y Q 



dx/dz~ 

8y/9z o 
9z/ 9z 



y 


z 



r 





+ 


(f-l)+l 0 0 

0 (f-l)+l o 

0 0 (f-1) 

+1 - 

r 

g o o 

o g o 

o o g 

L 6 J 


9x/ 9z 
9y/9z Q 



y y 



+ u 


y 


+ 


0 


9z/9z 

OJ 




0 


0 0 
f 0 
0 f 


3 x/ 3x^ 

9i/9y o 

9x/ 9z 

0 


X 

# 

X 

9y/ 9x Q 

9y / 3y Q 

9y/9z Q 

= 

y 

y 

9 z/ ax 

9z/ 9y Q 

9 * /9 V 


z 

# 

z 

— 


• # 

X 

(?o *0 g 

1 

(g-l)+l 0 0 

• » 

y 

-1- 




0 (g-i)+i 0 

z 

_ 


0 0 (g-i)+i 

— 


give the 36 partial derivatives of x, y> z, i, y, i, with respect 

X 0’ y 0’ Z Q’ *0* V V Finall y* t ^ ie formulas 
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9 Xq/ 9 c 

9 x q / 9 y 

9 x^/ az 

9 y Q / 9 x 

3y 0 /3y 

9 y Q / 3 z 

9 z / 9 x 

L 0 

9 z Q /ay 

9 z^/a z 

^9 x / 9 x 

9 x Q / 9 y 

9 x / 3 z 

3 y Q / 9 x 

3 y Q / 3 y 

3y Q / 9z 

_3 z^/ 3 x 

3 z Q / 9 y 

9 z / 9 z. 

8 Xq/ 3 x 

9 x / 9 y 

9 x / 3 z 

3 y Q / 9 x 

9 y 0 / 3 y 

9 y 0 / 9 z 

3 Zg/9x 

3 z Q / 9 y 

3 z /3 Z 

0 -j 

9 x q / 9 x 

9 x Q / 9 y 

3 x Q /az 

9 y Q / 9x 

3 y Q / 9 y 

9 y Q / 9 z 

9 z / 9 x 

9 z q / 9 y 

dz^/dz 


— rji 

ax/ax Q ax/ay o ax/az ' 


ay/9x Q 9y/ay Q 9y/9z Q 


9z/ax^ 9z/9y 0 9z/az 


9y/ 9x 


0 


9x/9x 
9y / 9x 


0 


0 


9z/9x 


0 


9x/ 0x 
9y/9x 


0 


8z/ 9x 


0 


0 


9x/ax ax/ay n ax/az 


9y/ 9y 


0 


0 


9 z/ 9 x^ 9 z/ 9 y 


0 


9y/ 9 z 


0 


9z/az 


0 


9i/3y. Q 

9 x/9z q 

3y/ 3y 0 

9y / 9 z q 

9z / 9 y.Q 

3z/ 3z 

U — 

9x/9y Q 

3x/az Q 

ay/ay 

3y/3 z 0 

3z/9y Q 

3z/3 z 

oJ 

V V 

V with 


• * 

respect to x, y, z, x, y, z. This concludes the computation of all 
quantities which are output by the subroutine. 
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4. 7 SOME PROPERTIES OF THE SUBROUTINE 


A Fortran 4 listing of the double-precision subroutine TWOBDY is 
given below. The compiled subroutine occupies less than 1500 single¬ 
precision storage locations on the IBM 7094. The minimum running 
time is a little over 10 milliseconds and occurs for the trivial case 
where the input value of ?/ is zero. The program will almost always 
run in less than 100 milliseconds for a non-zero value of sr where 
no input approximation for ip is available. A good input approximation 
for ip can reduce the running time to as low as 15 milliseconds. Ex¬ 
tremely unusual cases with no input approximation for ip can take as 
much as 1000 milliseconds, but the solution will always be obtained. 

Accuracy checks indicate that errors in the outputs are always 
close to a lower bound below which accuracy improvements are un¬ 
warranted. The lower bound occurs because of uncertainties in the 
inputs due to the floating-point representation of numbers. For example, 
the absolute error in x should not be reduced below about 


3 x 

9x o 

. |io" 

16 | 
x o> + 


ax 

9y o 

. |10" 

1—' 

O 

+ 

9x i| 
a, !■ > 

z o 

-16 

10 z 

0 

8x 

|. 10" 

16.1 | 
XqM 

a 

X 

I. 1 lo" 

16, 

4 

a x i 

1 10 *16. 1 

1 10 Z 0 I 

8i o 

a 

y o 

y o 

a z ' ’ 

0 


8 x i| t -l6 i 

■ 

|. | 10 Ml + 

X | 




10- 16 r 


because each of the eight inputs x , y^, z^, x^, y^, z^, ^ , and r is 
known to only about sixteen significant figures. Partial derivatives have 
been checked by multiplyiig the 6x6 inverse partial derivative matrix. The 
difference between the product and the unit matrix is then easily ex¬ 
plainable in terms of errors caused by the floating-point representation 
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of the partial derivatives. 

Results obtained in some limited applications of the subroutine 
have been very satisfactory. The author will sincerely appreciate 
any suggestions, criticisms, or comments from those using the 
Fortran program. A Fortran version for the IBM 360 should be 
available shortly. There are several programs available which 
solve the equivalent problem accurately and efficiently for particular 
cases of two-body motion. The main advantage of the program de¬ 
scribed here is that it offers equal accuracy and efficiency with com¬ 
plete generality as well. This can become extremely important in 
applications where a compact program is desired to handle many 
different types of two-body motion. In any case, the generality is 
provided without any sacrifice of accuracy, efficiency, or compactness. 
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ooooonooo 


4. 8 THE FORTRAN 4 SYMBOLIC LISTING OF THE SUBROUTINE TWOBDY 


SIBFTC TWOBDY 

SUBROUTINE TwOBDY(SOrTAU»MU»PSI»S»P»PI»PMU»P0MU»ACC»ACC0*R#R0) 

GENERAL SOLUTION OF TWO BODY PROBLEM WITH PARTIAL DERIVATIVES 
FORTRAN 4 DOUBLE PRECISION SUBROUTINE FOR IBM 7094 WITH IBSYS SYSTEM 
SEE APRIL 1965 ASTRONOMICAL JOURNAL FOR FORMULATION BY W. H. GOODYEAR 

CALLING SEQUENCE IS AS FOLLOWS 

CALL TWOBDY(SO»TAU»MUrPSI»S»P»PI»PMU»P0MU*ACC»ACC0»R»R0) 

DOUBLE PRECISION QUANTITIES IN CALLING SEQUENCE ARE AS FOLLOWS 
DOUBLE PRECISION SO(6)rTAU»MU»PSI 
l»S(6)»p(6»6)»PI(6r6> »PMU(6)»P0MU(6> »ACC(3)»ACC0(3)»R»R0 

C INPUTS 

C SO(1>#S0(2)»S0(3)=X0* YO * Z0=POSITION COMPONENTS AT REFERENCE TIME TO 
C S0(4)»S0(5)»SO(6)=XDO»YDO»ZDO=VELOCITY COMPONENTS AT REFERENCE TIME TO 
C TAU=TIME INTERVAL (T-TO) FROM REFERENCE TIME TO TO SOLUTION TIME T 
C MU=CONSTANT IN DIFFERENTIAL EQUATIONS (XDD>YoD»ZDD)=-MU*(X t Y»Z)/(R**3) 

C PSI=APPR0XIMATI0N FOR FINAL SOLUTION PSI OF KEPLER'S EQUATION 
C 

c OUTPUTS 

C PSI=GENERALIZED ECCENTRIC ANOMALY=SOLUTION OF KEPLERS EQUATION 
C S(l)fS(2)»S(3)=X»Y»Z=POSITlON COMPONENTS AT SOLUTION TIME T=T0+TAU 
C S(4)»S(5)*S(b)=Xu»YD»ZD=VELOCITY COMPONENTS AT SOLUTION TIME T=T0+TAU 
C P(I*J)=PARTIAL DERIVATIVE DS(I)/DS0(J) OF S<I) WITH RESPECT TO S0(J) 

C PI(I * J)=PARTIAL US0(I)/DS(J> WITH ROLES OF TO AND T REVERSED 

c pmui i impartial ds<i)/dmu of s(i> with respect to mu RLVtK5ED 

C P0MU(I)=PARTIAL USO(I)/DMU WITH ROLES OF TO AND T REVERSED 
C ACC<I)=-MU*S<I)/<R**3)=ACCELERATI0N COMPONENT AT SOLUTION TIME T 
C ACCO(I)=-MU*S0(I)/(R0**3)^ACCELERATION COMPONENT AT REFERENCE TIME TO 
C R-RADIUS AT TIME T-SQUARE ROOT OF(X**2+Y**2+Z**2) 

C R0=RADIUS AT TIME T0=SQUARE ROOT OF(X0**2+Y0**2+Z0**2) 
c 

C ADDITIONAL DOUBLE PRECISION QUANTITIES FOR COMPUTATION 

2»SIG0#ALPHA*PSIN.PSIP*A»APrC0rCl»C2»C3»C4>C5X3»SlfS2»S3#OTAU#DTAUN 
3»DTAUP#U»FMl»GfFD»GDM1 
C 

c start of initial computations 

C COMPUTE RADIUS RU=SQUARE ROOT OF(X0**2+Y0**2+Z0**2) 

S1=DMAX1<DABS(S0(1))»DABS(S0(2)) t DABS(SO(3))) 

S2=(S0(1)/S1)**2+(S0(2)/S1)**2+(S0(3)/S1)**2 

R0=2. 

10 R=R0 
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R0=(R+S2/R)*.f> 

IF(RO.LT.k) GO TO 10 
R0=R0*S1 

C COMPUTE OTHER PARAMETERS 

SIGO=SO(1)*SO<4)+SO<2>*SO(5>+SO(3)*SO<6) 
ALPHA=S0(4)**2+S0<5>**2+S0(6)**2-2.*MU/R0 
C INITIALIZE SERIES MOD COUNT M TO ZERO 
M=0 

C INITIALIZE BOUNDS PSIN AND PSIP FOR PSI OR SET PSI=0 IF TAU=0 
IF(TAU) 20#30»40 
20 PSIN=-l.U+38 
PSIP=0. 

DTAUN=PSIN 
DTAUP=-TAU 
GO TO 50 
30 PSI=0. 

GO TO 100 
40 PS1N=0. 

PSIP=+1.0+38 

dtaun=-tau 

DTAUP=PSIP 

C USE APPROXIMATION FOR PSI IF IT IS BETWEEN BOUNDS PSIN AND PSIP 
50 IF(PSI.GT.PSIN.AND.PSI.LT.PSIP) GO TO 100 
C TRY NEWTON’S METhOD FOR INITIAL PSI SET EQUAL TO ZERO 
PSI=TAU/R0 

C SET PSI=TAU IF NEWTON’S METHOD FAILS 

IF(PSI.LE.PSIN.OR.PSI.GE.PSIP) PSI=TAU 
C END OF INITIAL COMPUTATIONS 
C 

C BEGINNING OF LOOP FOR SOLVING KEPLER’S EQUATION 
C BEGINNING OF SER1LS SUMMATION 

C COMPUTE ARGUMENT A IN REDUCED SERIES OBTAINED BY FACTORING OUT PSI’S 
100 A=ALPHA*PSI*PSI 

IF(DABS(A)>Lt*1•) GO TO 120 

C SAVE A IN AP AND MOD A IF IT EXCEEDS UNITY IN MAGNITUDE 
AP=A 

110 M=M+1 
A=A*.25 

IF(DABS(A).G1»1») GO TO 110 
C SUM SERIES C5X3=3*S5/PSI**5 AND C4=S4/PSI**4 

120 C5X3=(l.+(l.+(l.+(l.+(l.+(l.+(l.+A/342.)*A/272.)*A/210.)*A/15b.) 
1 *A/110.)*A/72.)*A/42.)/40. 

C4 =(1.+<i.+(1.4<l.+<l.+<l.+(l.+A/306.)*A/240.)*A/182.)*A/132.> 
1 *A/90.)»A/56.)*A/30,)/24. 

C COMPUTE SERIES C3=S3/PSI**3»C2=S2/PSI**2»C1=S1/PSI»C0=S0 
C3=(.5+A*C5X^)/3. 
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C2= •b*A*C4 
Cl = l.+A*C3 
CQ= 1•4A*C2 
IF(M.LE.O) GO TO 140 

C DEMOD SERIES CO AND Cl IF NECESSARY WITH DOUBLE ANGLE FORMULAS 
130 C1=C1*C0 

C0=2.*CU*C0-1. 

M=M-1 

IF<M.GT.O) GO TO 130 

C DETERMINE C2»C3»C4.C5X3 FROM C0.C1.AP IF DEMOD REQUIRED 
C2=(C0-l,)/Ap 
C3=(C1-1.)/AP 
C4=(C2-.b)/Al- 
C5X3=(3.*C3-.b)/AP 

C COMPUTE SERIES Sl.S2.S3 FROM C1.C2.C3 
140 S1=C1*PSI 

S2=C2*PS1*PS1 
S3=C3*PSI*PSI*PS1 
C END OF SERIES SUMMATION 

C COMPUTE RESIDUAL DTAU AND SLOPE R FOR KEPLER’S EQUATION 
G=R0*Sl*SIG0*S2 
DTAU=<G+MU*S3>-TAU 
R=DABS(R0*C0+(SIG0*S1+MU*S2)) 

IF(DTAU) 200.300*210 
C RESET BOUND 
200 PSIN=PSI 
DTAUN=UTAU 
GO TO 22u 
210 PSIP=PSI 
DTAUP=DTAU 

C TRY NEWTON’S METnOD AND INITIALIZE SELECTOR N 
220 PSI=PSI-DTAU/R 
N=0 

C ACCEPT PSI IF IT IS BETWEEN BOUNDS PSIN AND PSIP 
230 IF(PSI.GT.PSIN.AND.PSI.LT.PSIP) GO TO 100 
C SELECT ALTERNATE METHOD OF COMPUTING PSI OR STOP ITERATIONS 
N=N+1 

GO TO (1.2.3.4.300).N 

C TRY INCREMENTING BOUND WITH DTAU NEAREST ZERO BY THE RATIO 4*DTAU/TAU 

1 IF(DABS(DTAUN).LT.DABS(DTAUP)) PSI=PSIN*(1.-(4.*DTAUN)/TAU) 
IF(DABS(dTAUP).LT.OABS(DTAUN)) PSI=PSIP*(1.-(4.*DTAUP)/TAU) 

GO TO 230 

C TRY DOUBLING BOUND CLOSEST TO ZERO 

2 IF(TAU,GT.O.) PSI=PSIN+PSIN 


94 



IF(TAU.LT.O.) psi=psip+psip 
GO TO 230 

C TRY INTERPOLATION BETWEEN BOUNDS 

3 PSI=PSIN+(PSIP-PSIN)*(-DTAUN/(DTAUP-0TAUN)) 

GO TO 230 

C TRY HALVING BETWEEN BOUNDS 

4 PSI=PSIN+(PSIP-PSIN>*.5 
GO TO 230 

C END OF LOOP FOR SOLVING KEPLER'S EQUATION 

C 

C COMPUTE REMAINING THREE OF FOUR FUNCTIONS FM1*GrFD>GDM1 
300 FM1=-MU*S2/R0 
FD=“MU*S1/R0/R 
GDMl=-MU*S2/R 

C COMPUTE COORDINATES AT SOLUTION TIME T=T0+TAU 
DO 310 1=1*3 

S(I)=S0(I> + (FM1*S0(I)+G*S0(I+3)) 

SU+3) = (FD*S0(I)+GDMl*S0(I+3) )+S0(I+3) 

C COMPUTE ACCELERATIONS 

ACC(I)=-MU*S(I)/R/R/R 
310 ACCO (I} =-MU*S0 (D/RO/RQ/RO 

C END OF COMPUTATION FOR COORDINATES AND ACCELERATIONS 

C 

C COMPUTATION OF PARTIAL DERIVATIVES 

C COMPUTE COEFFICIENTS FOR STATE PARTIALS 

U= S2*TAU+MU*(C4-C5X3)*PSI*PSI*PSI*PSI*PSI 

P(1»1)=-(FD*S1+FM1/R0)/RO 

P<1*2)=-FD*S2 

P(2»l)= FM1*S1/R0 

P(2*2>= FM1*S2 

P (1»3) = P(l»2) 

P(1*4)=-GDMl*S2 
P<2*3)= P(2»2) 

P(2»4)= G*S2 

P(3*1)=-FD*<C0/R0/R+1./R/R+1./R0/R0) 
P(3»2)=-(FD*S1*G0M1/R)/R 
P(4*1)=-P(1♦1) 

P(4»2)=-P(l»2) 

P(3»3)= P(3»2> 

P(3* 4)=-GDMl*Sl/R 
P(4*3)=-P<1*2) 

P(4r4)=-P(l»4) 

C COMPUTE COEFFICIENTS FOR MU PARTIALS 
P(l»5)=-Sl/R0/R 
P(2*5)= S2/R0 
P(3»5)= U/R0-S3 
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P(l»6)=-P(l,b) 

P<2»6)= S2/R 
P(3»6)=-U/R+S3 
DO 400 1=1,3 
C COMPUTE MU PARTIALS 

PMU <I)=-S <I> *P(2»5)+S <I+3)*P(3,5) 

PMU(I+3)= S<I>*P(l»5)+S<I+3)*P(2,5)+ACC<I)*P(3r5) 

P0MU< I)=-S0<i)*P(2,6)+S0U+3)*P<3»6) 

P0MU<I+3)= S0U)*P(1,6)+S0(I+3)*P(2»6)+ACC0(I)*P(3,6) 

C MATRIX ACCUMULATIONS FOR STATE PARTIALS 
DO 400 J=1,4 

PI(J»I)= P(J,D*S0(I)+P(j»2)*S0(I+3) 

400 PI(J»I+3)= P(J»3)*S0(I)+P(J»4)*S0(I+3) 

DO 410 1=1,3 
DO 420 J=1,3 

P(I»J) =S(I)*PI(1,J> +S(I+3)*PI(2,J) +U*S(I+3)*ACC0<J) 

P(I»J+3) =S(I)*PI<1»J+3)+S(1+3)*PI(2,J+3)-U*S(1+3)*S0(J+3) 

P(1+3,u) =S(I)*PI<3,J) +S(1+3)*PI(4,J) +U+ACC<I)+ACC0(J) 
420 P(1+3#J+3)=S(I)*PI(3,J+3)+S(1+3)*PI(4,J+3) -U+ACC (I) *S0(J+3) 
P(I»I) =P(I,I) +FM1+1. 

P(I»I+3) =P(I,I+3) +6 

P(I+3*I) =P(I+3#I) +FD 

410 P<1+3#I+3)=P(1+3#I+3)+GDMl+l. 

C TRANSPOSITIONS FoK INVERSE STATE PARTIALS 
DO 430 1=1,3 
DO 430 J=1» 3 
PI(J+3,I+3) = P(I,J) 

PI(J+3,I) =-P(I+3»J) 

PI<J,1+3) =-P(I,J+3) 

430 PI(J,I) = P(1+3,J+3) 

C END OF COMPUTATION FOR PARTIAL DERIVATIVES 
C 

C END OF PROGRAM - ALL OUTPUTS HAVE BEEN COMPUTEO 
RETURN 
END 
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